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Abstract 

This paper concerns pattern formation in a class of reaction-advection-diffusion 
systems modeling the population dynamics of two predators and one prey. We con¬ 
sider the biological situation that both predators forage along the population den¬ 
sity gradient of the preys which can defend themselves as a group. We prove the 
global existence and uniform boundedness of positive classical solutions for the fully 
parabolic system over a bounded domain with space dimension N = 1,2 and for the 
parabolic-parabolic-elliptic system over higher space dimensions. Linearized stabil¬ 
ity analysis shows that prey-taxis stabilizes the positive constant equilibrium if there 
is no group defense while it destabilizes the equilibrium otherwise. Then we obtain 
stationary and time-periodic nontrivial solutions of the system that bifurcate from 
the positive constant equilibrium. Moreover, the stability of these solutions is also 
analyzed in detail which provides a wave mode selection mechanism of nontrivial 
patterns for this strongly coupled system. Finally, we perform numerical simulations 
to illustrate and support our theoretical results. 
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1 Introduction 


One of the central problems in the study of ecological systems is to understand the spatial- 
temporal behaviors of population distributions of interacting species. Over the past few 
decades, many mathematical models have been proposed and developed to investigate 
the collective influence of individual’s dispersals on the spatial-temporal population dis¬ 
tribution at a group level. For example, reaction-diffusion equations have been applied 
to model the densities and spatial distributions of single or multiple interacting species in 
[13, 39,40, 4 ] etc. In particular, the formation of nontrivial patterns in these systems such 
as diffusion-driven heterogeneity or Turing’s pattern, especially those with concentration 
properties, can be used to describe the aggregation or segregation phenomena. 

In this paper, we study the following 3x3 reaction-advection-diffusion system 


Ut = V • (diVw — x u 4 > { w )'^ w ) + «i(l — u)u + fauw, 
v t — V • (d 2 Vn — £v<j>(w)Vw) + a 2 (l — v)v + favw, 

< U! t — d s Aw + 03(1 — w)w — fa\UW — /?32 vw, 

d n u = d n v = d„w = 0, 

k u{x, 0) = u 0 (x),v(x, 0) = v 0 (x),w(x, 0) = w 0 (x), 


x G f), t > 0, 

x G O, t > 0, 
x G t > 0, 
x G £2, 


( 1 . 1 ) 


where fl is a bounded domain in M ;V , N > 1 with smooth boundary dQ and unit outer nor¬ 
mal n. V is the gradient operator and A is the Laplace operator, u, v and w are functions 
of space-time variable (x, t ). di, otu i = 1,2,3, and fa, fai, i — 1,2, are positive constants. 
X and £ are also assumed to be positive constants and 4>{w) is a smooth function. 

(1.1) models the population dynamics of three interacting species subject to Lotka- 
Volterra kinetics, where u and v are population densities of two predators at location x 
and time t, and w is the population density of the prey. It is assumed that both predators 
consume the same prey species and disperse over the habitat by a combination of random 
diffusion and directed movement (prey-taxis) along the gradient of prey population den¬ 
sity, while the preys move in the habitat randomly. Diffusion rates di, i — 1,2,3, measure 
the intensity of random dispersals of the species. Here prey-taxis is the phenomena that 
predators with the ability to perceive the heterogeneity of prey distribution approach the 
patches with high preys density. The positive constants x and £ measure the intensity of 
the directed movement of each predator in response to the prey-taxis. The prey-dependent 
function 0 reflects the strength of prey-tactic movement of the predators with respect to 
the variation of prey population density. Various specific forms of 0 can be chosen, de¬ 
pending on the biological situation that one tries to model. In particular, xm0(w) > 0 
represents the biological situation that predator u forage the preys while x v M w ) < 0 
means that predator u retreat the habitat of the preys which can defend themselves as a 
group when the population density is high. See [ 5, 38, 64] and the references therein 
for more examples and further discussions about antipredation of preys through group 
defense. The population kinetics are of classical Lotka-Volterra type, where a* are the 
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intrinsic growth rates of the species which reproduce logistically, are the growth rates 
of the predators and are the death rates of the preys due to predation. 

In search of preys, prey-taxis is the process that predators move preferentially towards 
patches with high density of prey. All predators forage preys in accordance with spatial 
distributions of prey density. Though individual predator tends to forage the vicinity of 
recent captures, swarms of predators exhibit prey-taxis behaviors. This uneven searching 
effort may result in both larger predation rate and aggregation of predators where preys 
are abundant. For example, insects can find preys which are concentrated in some areas 
and sparse in others, though they lack long-distance sensory perception. See [ 16] and 
the references therein for detailed discussions and field experiments. This mechanism is 
the same as bacterial chemotaxis in which cellular organisms sense and move along the 
concentration gradient of stimulating chemicals in the environment [18, 19, 20, 28]. 

Various reaction-diffusion systems have been proposed to describe spatial predator- 
prey distributions under directed movements. In [ \6\, Kareiva and Odell proposed a mech¬ 
anistic approach, formulated as partial differential equations with spatially varying dis¬ 
persals and advection, to demonstrate and explain that area-restricted search does create 
predator aggregation. Since then, various reaction-diffusion systems have been proposed 
to model predator-prey dynamics with prey-taxis. Sapoukhina et al. [49] assumed that 
such directed movement is determined by the velocity variation (i.e. the acceleration). 
They investigated the effect of prey-taxis on the predator’s ability to maintain pest popu¬ 
lation below a certain economic threshold value. It is assumed in [ ,14, 56] etc. that the 

directed movement of predator is due to the advective velocity. We refer to [ 18, 28, 46, 57] 
etc. for the derivation or justification of (1.1). 

Reaction-diffusion models with prey-taxis subject to different population dynamics 
have been studied by various authors extensively. Ainseba et al. [2] established the exis¬ 
tence and uniqueness of weak solutions to prey-taxis models with volume-filling effect 
in prey-tactic sensitivity function. Global existence and boundedness of classical solu¬ 
tions are obtained in [16, 53], while nonconstant positive steady states are investigated in 
[58, 59, 61] and [ ] using Crandall-Rabinowitz bifurcation theory and fixed point the¬ 
ory, respectively. Lee et al. [ ] studied pattern formation in prey-taxis system under a 

variety of nonlinear functional responses. Travelling wave solutions for prey-taxis mod¬ 
els are studied by the same trio in [32]. Effects of prey-taxis on predator-prey models are 
investigated numerically in [ ] which suggest that both response functions and initial data 
play important roles in pattern formation; moreover, predator-prey models admit chaotic 
patterns when the prey-taxis coefficient is sufficiently large. We also want to mention that 
for predator-prey systems without prey-taxis, Okubo and Levin [ ] pointed out that an 

Allee effect in the functional response and a density-dependent death rate of the preda¬ 
tor are necessary for the pattern formation. See the book of Murray [ 0] for more works 
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on prey-taxis models. We also want to mention that there are also many works on the 
predator-prey models with density-dependent diffusion or the so-called cross-diffusions, 
such as [29, 30, 43, 48, 65] etc. 

All the aforementioned works are devoted to studying prey-taxis models with one 
predator and one prey. There are some works on two-predator and one-prey model with¬ 
out prey-taxis [22, 37, 45, 55]. Moreover, we refer the reader to [23, 51,5] etc. for the 
studies on the dynamics of competitive reaction-diffusion systems with cross-diffusion 
or advection. In [ ], Lin et al. considered (1.1) with % = £ = 0 over Q = (0, L) with 

or without diffusion. They investigated the global dynamics of all equilibria of the ODE 
system and obtained traveling wave solutions of the PDE system. In this paper, we con¬ 
sider system (1.1) with two predators and one prey, both predators foraging the preys 
prey-tactically. We are motivated to study the effect of prey-taxis and group defense on 
the spatial-temporal dynamics of (1.1), in particular, its positive solutions that exhibit 
stationary or time-oscillating spatial structures. 

There are several scientific goals of our paper and the rest part of this paper is orga¬ 
nized as follows. In Section 2, global existence of positive classical solutions are obtained 
for the full parabolic system in Theorem 2.6 over ID and 2D bounded domains and for 
the parabolic-parabolic-elliptic system over arbitrary space dimensions in Theorem 2.7 
respectively. We also prove that these classical solutions are uniformly bounded in time. 
In Section 3, we study the linearized stability of the positive equilibrium to (1.1). It is 
shown that prey-taxis destabilizes the positive equilibrium if there is group defense in 
preys and it stabilizes the equilibrium otherwise. Section 4 and Section 5 are devoted to 
the steady state and Hopf bifurcation analysis of (1.1) over Q = (0,L). Existence and 
stability of stationary and time-periodic spatial patterns are established in Theorem 4.2, 
Theorem 4.3 and Theorem 5.2. Though there have been many works devoted to the study 
of nonconstant steady states of the 2x2 prey-taxis systems, regular time-periodic spatial 
patterns are quite new to the literature. Finally, we perform extensive numerical studies in 
Section 6 to illustrate and support our theoretical findings. We would like to note that, C 
and Ci are assumed to be positive constants that may vary from line to line in the sequel. 


2 Existence and boundedness of global solutions 

In this section, we study the global existence and boundedness of classical positive solu¬ 
tions to (1.1). There are two well-established methods in proving the global boundedness 
for reaction-diffusion systems in the literature. One is to construct its time-monotone 
Lyapunov-functional and the other is to use Gagliardo-Nirenberg type estimate through 
Moser-ZT iteration. As we shall see in our mathematical analysis and numerical simula¬ 
tions, (1.1) admits time-periodic spatially inhomogeneous solutions for properly chosen 
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parameters and hence lacks time-monotone Lyapunov-functional and maximum princi¬ 
ple. Our proof of global existence and boundedness of classical positive solutions to (1.1) 
is based on the local theory of Amann [ , 6] and the Moser-Alikakos 17 iteration tech¬ 
nique [ ]. 


2.1 Local existence and preliminary results 

We first obtain the local existence and uniqueness of positive classical solutions to (1.1) 
and their extensibility criterion based on Amann’s theory [6]. To this end, we convert it 
into the following triangular form 



«i(l — u)u + (3\uw 
« 2(1 — v)v + 02 vw 
a 3 (l — w)w — /3 3 iuw — f3 32 vw 


( 2 . 1 ) 


with 


M 0 -xucj)(w)\ 

X>o = 0 d 2 - £v<t){w) , 

\0 0 d 3 J 

subject to nonnegative initial data m 0 , Vo, wq and homogeneous Neumann boundary con¬ 
ditions. Since all the eigenvalues of V 0 are positive, system (2.1) is normally parabolic, 
and we have from the standard parabolic maximum principles that u, v, w > 0 in Cl x M + . 
Moreover, the following results are evident from Theorem 7.3 and Theorem 9.3 of [ ], 
Theorem 5.2 in [6] and the standard parabolic regularity arguments. 


Theorem 2.1. Let Cl be a bounded domain in WL N , N > 1 with smooth boundary dCl. 
Assume that d l} oti for i = 1,2,3, and /),, f3 3i for i — 1,2 are positive constants. Sup¬ 
pose that the initial data (uq,vq,wq) G C 1 (f2) x C 1 (f2) x W 2,p for some p > N, and 
u 0 ,vo,w 0 > 0, ^ 0 in Q. Then there exist a constant T max e (0,oo] and a unique 
solution {u{x,t),v{x,t),w(x,t)) to (1.1) which is nonnegative on Cl x [0,T max ) such 
that (u(-,t),v(-,t),w(-,t)) G C([0, T max ), 77 1 (f2) x 77 1 (f2) x H 1 ( O)) and ( u,v,w) G 
C ioc ’ 2 (Cl x (0 ,T max )) 3 for any 0 < a < moreover, if sup sG(o t) || (u, v, w)(-, s)|| L °o 

is bounded for t G (0, T max ], then T max = oo, i.e., (u, v, w ) is a global solution to ( 1.1). 
Furthermore, ( u,v,w ) is a classical solution and ( u,v,w ) G C Q ((0, oo), C' 2( ^ 1_/3 ^(0) x 
C ,2 ( 1 -7)(f2) x (Cl)) for any 0 < a < /3 < 1. 


According to Theorem 2.1, the local solutions (u, v, w ) are global if their L°°-norms 
are bounded in time. To establish the global existence, we state some basic properties of 
the local solutions. 


Lemma 2.2. Suppose that all the conditions in Theorem 2.1 hold and («o, vq, wq) G 
C 1 (f2) x C' 1 (f2) x W 2,p for some p > N. Let ( u,v,w ) be the classical solution of (1.1) 
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over Cl x (0, T max ). Then there exist positive constants C\ and C 2 such that 


IIw(-, £)||l 1 ^) + ||u(-, ^)||z, 1 (n) < C\ — C(||mo||li, H^ollz, 1 , «i 5 &i, |^|), Vi G (0, T max ), 

( 2 . 2 ) 

and 

0 < w(x,t) <C 2 = C'(||w 0 ||oo,Q! 3 ),Vi G (0, T max ); (2.3) 

moreover, w(x, t ) G (0, 1) on Cl x (0, T max ) ifwo(x ) G [0,1] on (1 


Proof. We first prove (2.3). Let w(t) be the solution of 


dw(t ) 
dt 


a 3 (l — w(t))w(t),w( 0) = max vj 0 (x ), 

n 


(2.4) 


then it is easy to see that w(t) is uniformly bounded and it is a super-solution to the third 
equation of (1.1). Hence we have that w(x, t ) < w(t) from the maximum principle and 
this gives rise to (2.3). Moreover, if vj 0 (x) G (0,1), w(t) = 1 is a super-solution and this 
implies that w(x, t) G (0,1). 

We next prove the boundedness of ||u(-,f)||z,i and the same argument applies for 
IK^Ik. Integrating the first equation in (1.1) over Cl, we have from the boundedness 
of L oo that 

^ J u(x,t) < +/3i||u;(-,f)|| Il ooj J u(x,t)-a 1 J u 2 (x,t), 


here and in the rest of this section we skipped dx. In light of ( J Q u(x, t)) 2 < |0| f n u 2 (x, t), 
we have 


d 

dt 


u(x,t ) < («i + /3i\\w(-,t)\\ L oo) 




Solving this ordinary differential inequality gives 


2 


u(x, t) < max { f u 0 (x),\Cl\(l + — sup \\w(-, f)^^) ). 

! l in \ a i o<t<T max / J 


D 


In order to estimate L p -norms of u and v, we need to provide an a priori estimate on 
IIVHI L oo . The following lemma is due to the well-known smoothing properties of oper¬ 
ator — d 3 A + 1 and embeddings between the analytic semigroups generated by (e /A } />0 . 
We refer the reader to [21, 63] for references. 

Lemma 2.3. Suppose that Q is a bounded domain in R A , N > 1. Assume that all the 
conditions in Theorem 2.1 hold and let ( u,v,w) be the classical solution of (1.1) over 
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(0, T max ). Then there exists a constant C dependent on Htuollw 1 -®^) an< i -- such that for 
p e (1, oo) 

< c(l + sup (||«(-,s)|| L p + ||u(-,s)|| L p)Yvt G (0,T max ), (2.5) 

v se(o,t) ' 


for any q G (1, ifp G [1, N), q G (1, oo) ifp = N and q = oo if p > N. 

Proof We first write //’-equation into the following variation-of-constants formula 

wf,t) = e^ sA_1)t wo + [ e ( ' d3A ~ 1 ^ t ~ s ' ) h(uf,s),vf,s),w(- : s))ds, (2.6) 

Jo 

where h(u,v,w) = w + a 3 (l — w)w — f3 3 \uw — fJ 32 vw. Applying Lemma 1.3 in [ ] 

on (2.6) and using (2.3), we see that there exists a positive constant C such that for 1 < 

p, q < oo, 


\\w(-,t)\\wi«<c(l +e u{t s \t - s) s 2( p 9 ) (||m(-,s)|| L p + ||u(-,s)|| L p)/is^ , 

(2.7) 

where v is the first Neumann eigenvalue of —d 3 A + 1. On the other hand, we know from 
the gamma function that 


sup 

ts(0,oo) 


f e u< ' t s \t — s) 2 2 (p ^ds < oo, since ~ \ ~ ~ ~) > — 1, 

Jo 2 2 p q 


therefore (2.5) follows from (2.7). 


O 


In order to apply the Moser-iteration, it is sufficient to prove the boundedness of 
|| Vw(-, f)|| L 2 (jv+i) thanks to the quadratic decay kinetics in (1.1). 

Proposition 2.1. Let he a bounded domain in WL N , N > 1, and (u, v, w ) be the classical 
solutions of (1.1) over (0, T max ). If there exists a constant C\ such that 

||Vw(-, ^)IIl 2 ( jv + 1 ) < Cl,'it G (0, Tmax), 

then there exists a constant C 2 such that 


IK-,*)IU« + H^(‘, t)\\ l°° < c 2 ,it G (0,T max ). 

Therefore (u, v, w) is a global solution to (1.1) over 0 x (0, oo). 

Proof. For each p > 1, we have from straightforward calculations 

[ u p = [ u p ~ l V • (diVu — xucj)(w)Vw) + [ «i(l — u)u p + f>\u v w 


p dt 


in 


in 


' 1 d i (p - 1) 


pZ 


|VuS| 2 + 


2 X(P ~ 1) 


in 


P 


«20(u/)Vu 2 ■ Viu 


+ / «i(l — u)u p + /3iu p w. 


( 2 . 8 ) 






8 


Kc Wang, Qi Wang and Feng Yu 


To estimate the second integral in (2.8), we have from Young’s inequality that, for any 
constant e > 0, there exists C(e) >0 such that 

uiVu* ■ Vw < e|Vu 2 | 2 + e(M 2 )^ + C(e)\Vw\ 2(p+1) . 

On the other hand, for each p > 0 there exists C 0 > 0 such that a i (1 — u)u p + 3\v p w < 
— ^u p+1 + Co . For the simplicity of notations and without loss of our generality, we 
assume that HHh oo < 1 in light of (2.3). Choosing p = N in (2.8), thanks to the 
boundedness of ||Vw(-,£)|| L 2 (jv+i) we have 


1 d 

Ndt ' U 


N 


< 


< 


\S7u 2 |" + 


|Vm 2 | 2 + 


< 


4di(N — 1) 

N 2 

4di(iV-l) 

N~ 2 ... 

2 X (N-l)C(e) 

N 

Adi(N — 1) 2 X (N — l)e 


2 x(N - 1 ) 
N 

£,2 , *X(N-1) 


-Vw / m^ + CoIOI 
i z 


in 


N 


eu 


N+l 


+ e| Vu~ 


|Vtu| 


2(N+1) _ / N+l 


2 J 

| Vit 


< - 


N 2 ' N 

+ (*£-!)«_<*) r, +1+C[€) 

J 

‘Id \ (.V - 1) 


u 
n 

N l9 
2 Z 


+ ^( e ) 


N 2 


iv^r-^ 


JV+l 


+ C( e ); 


(2.9) 


where e is chosen to be sufficiently small. Then we can have from (2.9) that ||w(-, t) || L iv < 
C, Vf e (0, T max ) and this, in light of Lemma 2.3, implies that ||Vw(-, t)\\ L g < C for any 
q e [l,oo). 

In particular, we note that || Vu>(-, f)|| L 2 (jv+ 2 ) is bounded, therefore we can choose p = 
N + 1 in (2.8) and perform the same calculations as in (2.9) to show that ||w(-, t)\\ L N+i is 
uniformly bounded. Once again by applying Lemma 2.3, we can show that || Vw(-, t) ||l°° 
is bounded. 

Without losing the generality of our analysis, we assume that 

L°° < 1, and 

then we have from (2.8) 


/\ P< 


p dtt 


p z 


\Vui\ 2 + 2x{p ~ 1) [ u l|V«l|+ 0 ; f 
V Jn Jn 


u p — a\ 


U‘ 


,p +i 


where a\ = oil + jd\ sup t>0 || w(-, t) We can further estimate it 


12/V< 

pdt J n 


< 


4di(p - 1) 

p'2 


|VrC| 2 + 


2 xip ~ 1 ) 


P 


in 


e|Vu§| 2 + ^u p ) +a[ f u p 


4d 1 (p-l) 2x(p-l)e 

o ' 


< - 


2 di(p - 1) 

p2 


p 


( x(p- 1 ) 


|v “ ! n^v + “; 


w 




2di 


u p , 


n 


( 2 . 10 ) 
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by taking e smaller than —. Letting p go to oc and applying the standard Moser-Al ik akos 

PX 

-iteration [ ] on (2.10), we can show that there exists a constant C > 0 such that 
oo < C. 'it e (0,T max ). Similarly, we can prove the uniform boundedness of 
Z/°° • The proof of this proposition completes. □ 


2.2 A priori estimates of fully parabolic system for N < 2 

According to Lemma 2.3 and Proposition 2.1, in order to obtain the global existence of 
(1.1), it is sufficient to establish the boundedness of u and v in their //-norms for some 
p > N. Lemma 2.2 already gives L 1 boundedness from which global existence of (1.1) in 
ID follows. In this section, we restrict our attention to study (1.1) over 2D and we want to 
prove the boundedness of \\(u,v)(-,t)\\ L 2 . First of all, we introduce several entropy-type 
inequalities to estimate weighted functions ||wlnu|| L i, ||ulnu||x,i and || Vw||l 2 , etc. The a 
priori estimates rely on Gagliardo-Nirenberg interpolation inequality for which [ ,54] 

are good references. We now give the following lemma. 


Lemma 2.4. Let Vt C M 2 be a bounded domain with smooth boundary Oil. Suppose that 
all the conditions in Theorem 2.1 are satisfied and ( u , v. w ) are the local solutions to (1.1) 
on Q x (0, T max ). Then there exists a positive constant C independent of T max such that 

||wIn w||li(Q) + ||nlnu|| L i(n) + || V'iu|| L 2 (n) < C,Wt e (0 ,T max ). (2.11) 


Proof. Using u-equation of (1.1) and the boundedness of ||0(w)|| L °°, we have from inte¬ 
gration by parts and Young’s inequality 


— [ ulnu = I (Inu+ ].)«/ 

dtl o Jn J 


= / V ■ (di Vu — x u 4 > ( w ) Vw)(ln u + 1) + / (a\ — a\U + /9itu)-u(ln u + 1) 


in 


Vu 


= — (diVu — xuf(w)S7w) - -b / («i — ql\u + f3iw)u(\nu + 1) 

./o U Jo 


< — d\ 


|Vm| l 


-Wu 


+ X I ■ Vtu — ^ / u 2 \nu + Ci, 


’n 


u 


u 


< -- 

“ 2 


d\ f |Vu| 2 


’n 


u 


-f- / u 2 lnu + C 2 / \Wwf + C 3 , 


( 2 . 12 ) 


’n 


where Ci, C 2 and 6' 3 are positive constants. Similarly we obtain from the (—equation 


d 


d 2 f |Vn| 2 a 2 


— / v In v < -, 

dt J n 2 J n v 


^ / v 2 In v + C 4 / |Vw| 4 + C 5 . (2.13) 

4 In Jn 
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On the other hand, we have from straightforward calculations 

1 d 

2 dt 


[ |Vw | 2 = [ Vw-'Vwt — — [ Aww t 


in 


— / Aw(d 3 Aw + 0 ( 3(1 — w)w — fJ 3 \uw — /3 32 vw) 

Jn 

< —d 3 f \Aw \ 2 + I \Aw\(C e u + C 7 v + C$) 

Jn Jn 

A — 7 T [ \Aw \ 2 + Cg [ U 2 + C 10 [ V 2 + C u , 


(2.14) 


where CJ s are positive constants. In light of (2.3), we have from Gagliardo-Nirenberg 
interpolation that there exists C* > 0 such that 


||Vn;|| 4 L 4 <a(||Am ||| 2 + l). 

Multiplying (2.15) by gi and then adding it to (2.12) and (2.13), we have 

d 


(2.15) 


— ( / ulnu + v\nv + — |Vw| 2 ') 
dt V/o 2 / 


< _ <k 
~ 2 


|V«|- 


a 1 


u 


-2 / u 2 \nu + C 2 / |V W | 4 + C 3 


'O 


do 


_ l Vl f 

2 v 
d 3 gi 


a 2 


-A / n 2 lnu + C 4 / |VH 4 + C 5 


/n 


| Aiy | + C 9 / u + C 10 / v +C n 


< -- 

“ 2 


d\ f |Vu | 2 d 2 f |Vw | 2 


a 2 


— / tT In u - - ir In v 


, n u 

+ (C 2 + C 4 ) / |Vtu| 


n v 

4 d 3 g 1 


'O 


'O 


|Aw | 2 + C 12 


a 1 


a 2 


< - Ai / hi u \ v 2 Inv + (C 2 + C 4 ) / |Viu| 


d 3 g 1 


|Aw|“ + C 12 , 


where C) are positive constants. Choose /i 4 large such that > C*(C' 2 + 6' 4 + 1). We 
apply the fact that s In s < es 2 In s + C 0 Vs > 0, and Young’s inequality to obtain 

f ulnu+ f v\nv+^~ I |Vu’| 2 ') + f ulnu+ f vlnv + ^ f \\7w\ 2 < C, 
dt V Jq Jn 2 Jn 2 J a J u 2 J n 

(2.16) 

where e = min{^, ^}. Solving (2.16) through Gronwall’s lemma gives rise to (2.11). 

□ 

Lemma 2.5. Under the same conditions as in Lemma 2.4, there exists a constant C > 0 
such that 

IK-,*)IUa + IK - »*)IU a < cyt e (o,T max ). ( 2 . 17 ) 
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Proof. We have from the integration by parts, Holder’s inequality and (2.3) that 


1 d f 2 

2 dfJ n U 


= — di |Vu | 2 — x / “V ■ (u(p(w)'V'w) + / («! — apu + fd\w)u 2 

= — di f | Vu| 2 — ~ f u 2 ((f)(w)Aw + 0 / (u>)| Vw| 2 ) + f («i — + /3iw)u 2 

Jn 2 Jq Jo 

<~di f |V«| 2 + |(||M||i3||Aiy|| L 3 + ||«||i 3 |Vw| 2 i3 )-y f u 3 + C 1: (2.18) 


where to derive the inequality we have assumed that \\(p(w)\\ L oo, \\(J)'(w)\\ L °° < 1 without 
loss of our generality. 

To estimate (2.18), we apply Lemma 3.5 of [ ] with p — 3 to have that for any 

2 1 

e > 0, there exists a constant C(e) > 0 such that ||m|| L 3 < e|| Vm||£ 2 \\u In u\\ 3 Ll + 

1 

C(e)(||wln-u|| L i + llwll^i). This fact and (2.11) imply that 

IMIls < (e||V u|| 2 2 + C 2 ) C (2.19) 


Thanks to the Gagliardo-Nirenberg inequality [5 Z ] and (2.11), we have 

||Aw|| L 3 < C 3 (|| VAw||f 2 || Vu;||| 2 + ||VHW < C 4 (||VA W ||| 2 + l), (2.20) 


and 



llvHlie < CsfllvAHlilMll- + |Mli~) < c 6 (||va w ||! 2 +1). 

( 2 . 21 ) 


Combining (2.19) with (2.20) and (2.21), we have from Young’s inequality that 

\\u\\ 2 l3 \\Aw\\ l3 < C 7 (e||V«||i 2 + Co) § (||VA W ||| 2 + l) 
<C 8 (6||V«||i 2 +C 0 )(||VAu;||l 2 + l) 
<| l I|v«H! 2 + 6||va w ||| 2 + c 9 


and that 


Mli. 


|Vw| 5 


< y II Vm|| 2 2 + e||VAw||^ 2 + Ci 0 , 
L3 Z X 


therefore (2.18) gives us 


1 d 

2 dt 


u 2 + 


di 

2 


|Vm| 2 < 2e / |VA w\‘ 


in 


OL\ 

t 



+ Cn. 


( 2 . 22 ) 














12 


Kc Wang, Qi Wang and Feng Yu 


By the same arguments, we can show that 

| Vn | 2 <2 el | VAie| 


1 d 

2 dt 


v 2 + 


d 2 


2 a 2 


V 3 + C\ 2 - 


(2.23) 


in 


On the other hand, we operate V to w-equation in (1.1). Then it follows from Young’s 
inequality and integration by parts that 

I \Aw\ 2 — I AwAw t = — f VA w-X7w t 

2 dt Jq J n 

= — VA w ■ V (d^Aw + 03(1 — w)w — /? 3 i uw — ferw) 

Jn 

= — dz / |VAu>| J + / \7Aw ■ (a^Vw — 2a 3 u?Vu> — ^\uVw) 


VA w ■ (—/ 3 3 iw'Vu — ( 3 3 2 fVru — feuiVw) 


d 3 


<-y / |VAu >| 2 — a 3 / |Au >| 2 + C 13 / \Vu\ 2 + C- 


14 


|VU | 2 + C; 


15 


|vh ! 


’n 


+ Oi 6 / m 2 |Vu7| 2 + Ci 7 / n 2 |Vuf 


d 


<(—V + 2e) / |VAiu| 2 — a 3 / |Atu| 2 + 2C' 13 / |Vm | 2 + 2C' 14 / |Vif + C 18 , 


(2.24) 


where we used the inequalities 


Ci 6 / u 2 |Vre | 2 < C 16 ||«|| 2 3 |Vw | 2 < C7 13 ||V«||i 2 + e||VA ^|| 2 2 + C7 19 

./O i 3 


and 


C 17 / n 2 | Vuf < C ' 16 ||^||| 3 | Vuf < C' 14 ||Vn || 2 2 + e||VA ^||| 2 + C 20 . 

Jn L3 

Multiplying (2.24) by sufficiently small with e = we add it with (2.22) and 

(2.23) to have 

d (l 


dt (2 L “ 2 + 2 L v2 + f L 1 Am|2 ) + ( L “ 2 + !/ + “ m L 1 Am,|2) - c - 

(2.25) 

Therefore (2.17) follows from (2.25) thanks to Gronwall’s inequality. □ 


2.3 Existence and boundedness of global solutions 

We now present our main results on the global existence and uniform boundedness of 
positive classical solutions to ( 1 . 1 ). 
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Theorem 2.6. Let Pi be a bounded domain in M. N , N < 2. Suppose that all the pa¬ 
rameters in (1.1) are positive and (j) E C 3 (M: R). Then for any nonnegative initial data 
(uq, Vo, wo) E C 1 (f2) x C 1 (f2) x W 2 ,p (Pl), p > N, there exists a unique triple (u, v, w) of 
nonnegative boundedfuncstions in C°(f2 x [0, oo)) fl C 2,1 (f2 x (0, oo)) which solve (1.1) 
classically in 13 x (0, oo). 

Proof. According to Theorem 2.1 and Lemma 2.3, we only need to show that ||(w(-,t),n(-,t))|| 
is uniformly bounded for all t G (0, T max ), then T max = oo and the existence part of The¬ 
orem 2.6 follows. Moreover, one can apply the standard parabolic boundary //-estimates 
and Schauder estimates in [ ] to verify that u t ,v t , w t and all spatial partial derivatives of 

u, v and w up to second order are bounded on PI x [0, oo), therefore (u, v, w) have the 
regularities stated in Theorem 2.6 

Thanks to Lemma 2.3 and Lemma 2.5, we have that ||Vw(-, t )\\lp is bounded for any 
p G (1, oo), then the uniform boundedness of ll(u,u)M)ll L oo follows from Proposition 
2.1 . Together with (2.3), we conclude the proof of Theorem 2.6. □ 


Our proof of the boundedness of ||(u, n)||i<x» hence the global existence of (1.1) is 
restricted to the domain with dimension N — 1,2 due to technical reasons. It is well 
known that, in order to prove the boundedness of ||Vtu||z,oo, it is sufficient to prove the 
//-boundedness of u and v for some p > 2(N + 1) for a wide class of reaction-diffusion 
system with biological relevance, for instance the Keller-Segel chemotaxis models with 
no cellular growth. The logistic decay in (1.1) contributes so one only needs to prove 
//-boundedness of u and v for some p > N. Whether or not the logistic dynamics 
are sufficient to prevent finite or infinite time blowups for (1.1) over higher dimensional 
domains is unclear in the literature and it is beyond the scope of this paper. However, we 
can show the global existence and boundedness for the following parabolic-parabolic- 
elliptic system of (1.1) in the following system regardless of space dimensions 


u t — V • (d]Vw — x u< f > ( w )'^ w ) + «i(l — u)u + (3iuw, 
v t — V ■ (d 2 Vv — £v<f>(w)Vw) + a 2 (l — v)v + /3 2 mu, 

0 = d 3 Aw + « 3(1 — w)w — /?3i UW — /?32 vw, 
d n u = d n v = d n w = 0 , 

k u(x, 0) = u 0 (x),v(x, 0) = v 0 (x),w(x, 0) = w Q (x), 


x E Pi, t > 0, 
x G Pi, t > 0, 
x E Pi, t > 0, 
x G DPI, t > 0, 
x E Pi, 


(2.26) 


To be precise, we have the following results. 


Theorem 2.7. Let Pi be a bounded domain in M /V , N > 3 with smooth boundary DPI. 
Suppose that (j>'(w) > 0, Ww > 0. Then (2.26) has a unique classical solution ( u,v,w) 
which exists globally in time and satisfies the following estimate with a positive constant 
C 
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Proof. By the same analysis as in the proof of Proposition 2.1 , we shall only need to show 
the uniform boundedness of ll«M)IU oo and ||r>(-, f)|| L oo. Similar to (2.8), we can obtain 
from the integration by parts 


ld [ p 

—— / vF 

P at J n 
_ 4di(p - 1) 

p2 

_ 4 di(p - 1 ) 

p2 


2,2 , XiP~ !) 


|VW5| 2 + 

\Vu^\ 2 - 


p 


x(p ~ 1 ) 


P 


X7u p ■ ((fi(w)'Vw) + / («i — aiu — f3\w)u p 


u p [cj)'{w) |Vru| J + 4>(vj)A w 


+ / («i — oiiu — /3iw)u p 


< 


4rfi(p - 1) 

p2 


| VvS 


x(p - 1 ) 

p 


u p cj)(w)Aw + / («i — ot\u — /3iw)u p . 




(2.27) 


In light of the pointwise identity d 3 Aw + a 3 (l — w)w — ffnuw — P^vw = 0, we have 


Id f . 

pdt U 


in 


< 


4di(p-l) /■ |V7 ^2 |2 x(p~ 1) 

d 3 p 


P Jn 


Vu 2 


(2.28) 

u p (p(w)(a 3 w(w — 1) — /3 3 ±uw — /3 3 2vw) 


+ / («i — aiM — f3iw)u p 


< 


4di(p - 1) 


P 7o 


|V«2| 2 -Ci / m p+1 + C 2 , 


(2.29) 


where Ci and C 2 are positive constants that depend on the system parameters and ||u>||l°°- 
Applying the Moser-iteration gives rise to the boundedness of IK->f)ll LjOO . By the same 
way we can prove the boundedness of ||u(-, The proof of this Theorem completes. 

D. 


The assumption (j>'{w) > 0 corresponds to the biological situation that the intensity of 
the directed dispersals of predator species increases as prey density increases, therefore 
there is no group defense in prey species. We have to leave it for future works on the 
global existence of (2.26) and its fully parabolic counter-part over higher-dimensions. 


3 Linearized stability of homogeneous equilibrium 

From the viewpoint of mathematical modeling, it is interesting and meaningful to investi¬ 
gate (1.1) for its nontrivial solutions describing the spatial distributions of the interacting 
species over the habitat. We shall see from the mathematical analysis that the formation 
of such nontrivial patterns is due to the joint effect of prey-taxis and sensitivity function. 
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To illustrate the effect of prey-taxis on the pattern formation and for the simplicity 
of calculations, we restrict our attention to (1.1) over one dimension and consider the 
following system 

{ u t = {dyu! — xu4>(w)w'y + «i(l — u)u + j3\uw 1 x € (0, L), t > 0, 
v t = (d 2 v'— £v(f)(w)w'y + ot 2 {l - v)v +/3 2 vw, xe(0,L),t>0, 

w t = d 3 w" + 0 : 3(1 — w)w — /? 3 i uw — (3 32 vw, x G (0, L), t > 0, 
u'(x , t) = v'(x, t ) = w'(x , t) = 0, x — 0, L, t > 0, 


where ' denotes the derivative taken with respect to x and all the parameters in (3.1) are 
the same as those in (1.1). 

We can find that (3.1) has six equilibria and one of them is 


(u,v,w) 



P1 - * @2 - _\ 

— w, 1 H- w, w ) 

(y,\ ot 9 / 


ft. 

-T" 

O 2 


a 3 ~~ /?31 — /?32 

+ M3i + Ma’ 

0 Q!l Q2 


which is positive if and only if o 3 > p 3 1 + /) 32 . To explore the existence and stability of 
stationary and oscillatory nonconstant positive solutions to (3.1), our starting point is the 
linear stability of (u, v. w) and we shall assume that it is positive from now on. 

Linearizing (3.1) around the constant equilibrium (u,v,w) and letting ( u,v,w ) = 
(u, v, w) + (U, V, W), U , V, W being small perturbations from (u, v , w), we obtain the 
following system of ( U , V, W) 


r U t tt (dif/' - xu(l>{w)W'y - ayaU + fauW, xe (0, L),t > 0, 

V t » (d 2 V - iv<j)iw)W'y - a 2 vV + fovW, xe(0,L),t>0, 

W t ~ d 3 W" — a 3 wW — /3 3 iwU — f3 32 wV, x e (0, L), t > 0, 

k £7 , (x)=y , ( a; ) = ^ , ( a; )=0» z = 0,L,f>0. 

According to the standard linearized stability principle ([5 ] e.g.), the stability of (u. v, w) 
is determined by eigenvalues to the following matrix 

/-^i(x) 2 _ a iu 0 x«0(w)(^) 2 + Piu\ 

0 -d^tff-a 2 v ^v^w)^) 2 + f3 2 v . (3.2) 

\ -P31W -P32W ~ 4 (x ) 2 - “ 3 ^ / 

We have the following result. 


Proposition 3.1. Assume that all the parameters in (3.1) are positive and the condition 
q,3 > /3 3 i + /3 32 is satisfied. Suppose that oiw) < 0 , then the positive equilibrium (u, v. w) 
is locally asymptotically stable if x < Xo and it is unstable ifx> Xo> where 


Xo = min {xt , Xk } 

fce N+ 


with 


PwvHi , _ H 1 H 2 H 3 + H 2 (3 3 i(3iuw + H 1 f3 32 f3 2 vw 
@ 3 iuH 2 H 2 /3 31 uw(j)(w)(!f ) 2 


(3.3) 

(3.4) 
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and 


H _ H 2 H 2 + H X H 2 + H 2 H 3 + II X II 2 , + Hi 11, + H 2 H 2 + 2// i // 2 // 3 

Xk 


(H x + + (H 2 + H 3 )f 32 f 2 v (H 2 + H 3 )p 32 v 


(- H 1 + H 3 )p 31 uci ) (w)(!fy 


(H x + H 3 )fd 31 u 


e, 


where 


TT 1 ( \ 9 — TT 7 / \ 9 — TT 7 / ^77 . O _ 

H\ = di( —) + «iM, Ho = d 2 (—) + at 2 v, H 3 — d 3 (—) +a 3 w. 


(3.5) 


Proof. The characteristic equation for an eigenvalue a of the stability matrix (3.2) is 


o -3 + mix, k )° 2 + Vi(x, k)a + mix, k ) = 0, 

where 

mix, k) =(di + d 2 + d3)(-y-) 2 + Oi\u + a 2 u + > 0, 

±j 

m(x,k) = ( rf i(^ r ) 2 + “ 1 '“) (^ 2 (^-) 2 + 02 ^ + [dii^-f + aiti'j (d 3 {^-) 2 + a 3 u>) 

+ (dni^) 2 + a 2 u) (dzi^) 2 + cy 3 w^j + (/?3l/3lW + p32l32v)w 
+ (folUX + /?32 V^wfiw)^) 2 , 

and 



/3 3 iuwcl)(w)(^-) 2 x + (dii^-) 2 + a!U^f3 32 vw(j)(w)(^) 2 £. 


By the principle of the linearized stability (Theorem 5.2 in [52] e.g.), (u, v, w) is asymp¬ 
totically stable with respect to (3.1) if and only if all eigenvalues of the matrix (3.2) have 
negative real part, then according to the Routh-Hurwitz conditions, or Corollary 2.2 in 
[16], the constant equilibrium (u,v,w) is locally asymptotically stable with respect to 
(1.1) if and only if the following conditions hold for each k <G N + 


mix, k) > 0, rh{x, k) > 0, and mix, k)m{x, k ) - mix, k) > 0, 

while («, v, w ) is unstable if one of the conditions above fails for some k e N + . Note that 
we always have that mix, k) > 0 for each k e N + ; moreover, //, (y, k) > 0 if mix, k) > 0 
and mix, k)mix, k) — mix, k) > 0, therefore (w, v, w) is unstable if there exists k e N + 
such that either mix, k) < 0 or //, (y, k)r/ 2 (x, k ) — mix, k ) < 0. On the other hand, since 
4>iw) < 0, it follows from straightforward calculations that 

Voix, Ac) < 0 if and only if x > Xk, 
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and 

7i(x, k)r} 2 (x, k) - Tj 0 (y, Ac) < 0 if and only if y > Xk, 

therefore the constant solution (u, v. w) is unstable if there exists k G N + such that X > 
Xk or X > Xk, i.e., if y is larger than the minimum of yf an( l Xk over k G N. Similarly 
we can show that ( u,v,w) is locally asymptotically stable if y < y 0 . This finishes the 
proof. □ 

According to Proposition 3.1, (u, v, w) becomes unstable as y surpasses the threshold 
value yo if < 0 which we shall assume from now on. We would like to point out that 
biologically 4>(w) < 0 describes situation that a huge amount of preys can aggregate to 
form group defense and keep predators away from the habitat when the prey population 
density surpasses w. This amounts to a switch from prey-attraction to prey-repulsion in 
the predation. See [ , 15, 38, ( ] and the references therein for the detailed description of 
prey group defense. 

We will see in our coming mathematical analysis and numerical simulations that 
(u, v, w) loses its stability to time-periodic patterns through Hopf bifurcation if y = Xk 
and to stationary patterns through steady state bifurcation if y = y f. Therefore we have 
used the indices H and S to denote Hopf and steady state bifurcation respectively. Here 
and in the rest part of this paper, we study the effect of prey-taxis on the formation of 
nontrivial patterns to (3.1). Without losing the generality of our analysis, we treat y as the 
variable parameter and fix all the rest parameters, while similarly, we have that Propo¬ 
sition 3.1 also holds for £ > £ 0 = min feeN +{£jf,£f}, where £jf and £f are functions of 
X- 

Corollary 1. Suppose that oiw) > 0 and all the rest conditions in Proposition 3.1 are 
satisfied, the homogeneous equilibrium (u, v, w) is locally asymptotically stable ifx > Xo 
and it is unstable ifx< Xo> where 


Xo = max{yf,yf}. 
fee N 

When 4>(w) > 0, both yf and y.k arc negative for any k G N + hence y 0 < 0, 
therefore (■ u,v,w ) is always stable according to Corollary 1 since y > 0. This result 
corresponds to the widely held belief ([ ] e.g.) that prey-taxis stabilizes homogeneous 

equilibrium and inhibits the formation of spatial patterns for one-predator and one-prey 
system. It states that the same holds true for two-predator and one-prey model as long as 
4>(iu) > 0. However, if o{vj) < 0, the prey-taxis destabilizes homogeneous equilibrium 
which becomes unstable as y surpasses yo - Therefore, in order to investigate the formation 
of nontrivial patterns in (3.1), we shall assume that 4>(iu) < 0 in the coming analysis. 
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Remark 3.1. As we shall see in our coming analysis, the occurrence of Hopf or steady 
state bifurcation at (u, v, w) depends on whether xo in (3.3) is achieved at min fceN + x£ or 
min feeN + Xk- We divide our discussions into the following cases: 


(1). If Xo = Xk 0 < min fceN+ Xk’ then Vo{Xo, fa) = 0 and the eigenvalues of (3.2) at xo 
are 


<rf{Xo,ko) = 0,of j3 (xo,fco) 


-m(xo, fa) ± VvKxo, fa) - 4yi(xo, fa) 

2 


Since (u, v, w) is unstable for all x > Xo, we know that (3.2) has at least one eigen¬ 
value with positive real part when x = Xk - k fa k 0 . This fact will be applied in the 
proof of the stability of steady state bifurcating solutions to (3.1) around (u, v, w, x£ )• 
In particular, it implies that the only stable bifurcating solutions must be on the 
branch around ( u,v,w,Xk 0 ) that turns to the right while all the rest branches are 
always unstable. Moreover, Hopf bifurcation does not occur at (u, v, w) in this case. 


(2). Ifxo = xZ < min fceN+ x£, then rjo(Xo, h) = rji(xo, ki)rfr(xo, h) and the eigenval¬ 
ues of (3.2) at xo are 

o-f (xo, fci) = -?? 2 (xo, fa) < 0, af 3 (x 0 , fa) = ±y/-‘n 1 (xo,fa). 


Similarly we can show that the steady state bifurcating solutions around (u,v, 
w, Xk) are always unstable for all k G N + . Moreover we see that rp (xo, fa) > 0 and 
(3.2) has three eigenvalues: erf (xo, fa) = -Viixo, fa), a^ 3 (xo, fa) = F^pfaxo, fa)i. 
This indicates the possibility of a Hopf bifurcation and the emergence of time-periodic 
spatial patterns in (3.1) when x = Xk, ■ To prove this claim, we argue by contradic¬ 
tion and assume that pi(xo, fa) < 0, therefore cr^ix o> ^i) — y—ViiXo, fa) > 0 and 
Re(a- 2 ) > 0 if x is sightly smaller than Xk,, however, this indicates that ( u,v,w ) is 
unstable for x < Xk, — Xo. which is a contradiction. Our numerical simulations 
in Section 6 support the existence of Hopf bifurcation and time-periodic patterns to 
(3.1). 


(3). If Xo — Xk (l — Xkc T'l this case, (3.2) has three eigenvalues ai = — r/ 2 (x> k) < 
0, a -2 = cr 3 = 0, and linear stability of the (?./, v, w) is lost since there are two zero 
eigenvalues. This inhibits the application of our steady state and Hopf bifurcation 
analysis which requires the null space of (3.2) to be one-dimensional. Therefore we 
assume that xf fa Xk ■ VA: G N + ' in our bifurcation analysis. 

In general, it is not obvious to determine when case (1) or case (2) in Remark 3.1 
occurs. However, if the interval length L is sufficiently small, we have 

v s _ d\d^ /kny did 3 + d 2 (d 1 + d 3 ) + d\ rki r\ 2 _ H + 

Xk ~ p 31 uwfiw) \ L ) /3 3 i uwfiw) V L ) ~ Xk ’ ’ 
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since 4>(w) < 0. This implies that for small intervals, xo = min fceN + xf = xf and 
(u,v,w) loses its stability to the steady state bifurcating solution as x surpasses Xo- 
Since k 0 = 1, the bifurcating solution has a stable wave mode cos ^ which is mono¬ 
tone in x. Moreover, we shall observe from our numerics that k 0 is increasing or non¬ 
decreasing in L. Therefore, small domain only supports monotone stable solutions, while 
large interval supports non-monotone solutions, at least when x is around xo • Actually, if 
(u(x), v(x), w(x)) is an increasing solution to (3.1), (u(L — x),v(L — x),w(L — x)) is a 
decreasing solution, then one can construct non-monotone solutions to (3.1) over (0, 2 L), 
(0, 3/.),...by reflecting and periodically extending the monotone ones at x = L, 2 L, 3 L,... 

4 Nonconstant positive steady states 

This section is devoted to studying nonconstant positive steady states to system (3.1), i.e., 
nonconstant solutions to the following system 

{ (diu' — xu ( j)(w)w , y + «i(l — u)u + /3iuw = 0, x G (0, L ), 

{d 2 v' - £v(j)(w)w')' + Q! 2 (l - v)v + /3 2 vw — 0, xe (0,L), 
d^w" + 03(1 — w)w — /? 3 i uw — flavin = 0, x e (0, L), 
u'(x) = v'(x) = w'(x ) =0, x = 0, L, 

where u, v and w are functions of x and all the parameters are the same as those in (3.1). 
We assume that a 3 > /3 3 i + /3 32 so that ( u , v , w) is the unique positive equilibrium to (4.1). 

In order to look for nonconstant positive solutions to (4.1), we shall perform steady 
state bifurcation analysis at (u, v, w). When 0 ( 10 ) < 0, we already know that prey-taxis 
X destabilizes ( u,v,w ) which becomes unstable when % surpasses xo, therefore we are 
concerned with the conditions under which the spatially inhomogeneous solutions emerge 
through bifurcation as x increases. We refer these as prey-taxis induced patterns in anal¬ 
ogy to Turing’s instability. 

4.1 Steady state bifurcation 

To apply the bifurcation theory of Crandall-Rabinowitz [9, 10] with x being the bifurca¬ 
tion parameter, we introduce the spaces 

X = {we H 2 ( 0, L)\w\0) = w'(L) = 0}, y = L 2 ( 0, L) 

and convert (4.1) into the following abstract equation 

J-(u, v, w, x) = 0, ( u , v, w, x) £ X x X x X x M, 


where 


F(u,v,w,x) 


(diu' — xu(t>{w)w')' + «i(l — u)u + f3iuw 
(d 2 v' — £v(j)(w)w , y + a 2 (l — v)v + / 3 2 vw 
d 3 w" + 03(1 — w)w — / 33 iUW — flavin 


(4.2) 
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It is easy to see that F{u, v, w, x) — 0 for any x € K and T : A’xA’xA’xl-^^x^xJ 
is analytic. Moreover, for any fixed ( u,v,w ) E X x X x X, the Frechet derivative of T 
is given by 

DJ r (u,v,w,x)(u,v,w) 

( dpi!' — x(ucj)(w)w' + uw(j)'{w)w' + u(j>{w)w')' + DT\ \ 

d 2 v" — Sfv<f(w)w' X vw<f'{w)w' X v<f(w)w')' X D!F 2 J, (4.3) 
d 3 w" - (3 3 iwu - /3 32 wv + (a 3 - 2 a 3 w - f3 31 u - P 32 v)w ) 

where DT\ = (an — 2a\u + f3iw)u + [3\uui and DT 2 = (a 2 — 2a 2 v + /3 2 ui)v + /3 2 fw. 

We collect the following facts about X. 

Lemma 4.1. DJ r (u,v,w,x){ u -i v - l w ) 'XxXxXxW—^yxyxy is a Fredholm 
operator with zero index. 

Proof. We denote u = (u, v, w) T and rewrite (4.3) as 

DJ r (u,v,w,x)( u , v , w ) = -4o(u)u" + F 0 (a:,u, u), 


where 

fdi 0 -xuc/)(w)\ 

Ao — 0 d 2 -£vcj)(w) 

\ 0 0 d 3 ) 

and 


( —x{u(j){w)w' + uwf (w)w')' + (■ u4>(w))'w' + DT A 
*-£(v<j>(w)w' + vw4>'{w)w')' + ( vf{w))'w' + DX 2 J , 

~/3 31 wu - (3 32 wv + (a 3 - 2a 3 u> - f3 31 u - (3 32 v)w J 

therefore operator (4.3) is elliptic since all eigenvalues of Aq are positive. According to 
Remark 2.5 (case 2) in Shi and Wang [50] with N — 1, Ao satisfies the Agmon’s condition 
(see Theorem 4.4 of [ ] and Definition 2.4 in [50]). Therefore, Dy{u, v, w, x)(w, v, w ) is 
Fredholm with zero index due to Theorem 3.3 and Remark 3.4 of [ ]. □ 


To seek non-trivial solutions of (4.1) that bifurcate from equilibrium (u, v, w), we first 
check the following necessary condition 

Af(DX(u,v,w,x))f{0}, (4.4) 

where A f denotes the null space. Taking (u, v , w) = (u, v, w) in (4.3), we see that the null 
space in (4.4) consists of solutions to the following problem 

" d j u" — xuf(w)w" — ci\uu + /3iuw = 0, x E (0, L), 
d 2 v" — £v(f)(w)w" — a 2 vv + / 3 2 vw = 0, x E (0, L), 
d 3 w" — f3 3 iwu — fd 32 wv — a 3 ww = 0, x E (0, L), 

^u'(x) = v'(x) = w'{x) = 0, x = 0,L. 


(4.5) 
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In order to verify (4.4), we substitute the following eigen-expansions into (4.5) 


OO y OO y OO 7 

, KTTX K7nE /C7TX 

= ^tkCOS—j—TVix) = 2_^s k cos——,w(x) = 2 JkCOS——. 


k =0 


k=0 


k =0 


t k , s k and r k constants and collect 


0 

' kn\ 2 



yu^w)^) 2 + /3iu\ /4 N 

— \ ( k7T \ 2 


-d 2 ( f f) 2 - a 2 v Zvcfiiw)^) 2 + /3 2 v I I s/c ) = I 0 


(4.6) 


~/3 32 w 


-d 3 (^) 2 - a 3 w 




k = 0 can be easily ruled out since a 3 > /3 31 + /? 32 . For each k G N + , (4.5) has nonzero 
solutions (4, Sfc, r/,.) if and only if the coefficient matrix of (4.6) is singular or equivalently 


s /? 3 2 vHi H\H 2 H 3 + H 2 f5^\f3\uw + Hif3 32 /3 2 vw 

X = Xk = ~ o _-TT Z ~ 


fd 3 iuH 2 


-4 2 /3 3 iMtu0(w)(¥) 2 


(4.7) 


where Hi, IL> and // 3 are given in Proposition 3.1. Note that x k i n (4-7) is the same as 
(3.4). x k > 0 if and only if 


£< 


H 1 H 2 H 3 + H 2 /3 31 /3iuw + Hi/3 32 f3 2 vw 
Hi(3 32 vw<P(w)(f) 2 


Condition (4.4) is satisfied if % = x.f and J\f(DJ r (u, v, w, x k )) — span {(u k ,v k ,w k )} 
which is of one dimension 


„ knx _ knx _ knx 

u k = P k cos —j~, v k = Q k cos —j-i w k = cos —j—, 


where 


P k = 


As 2 £W>M(^) 2 + /3 2 u d 3 (^f-) 2 + a 3 w 


An d 2 (>fy + a 2 v 


An w 


(4.8) 


(4.9) 


and 


^ k d 2 ( k jy + a 2 v ' 


(4.10) 


Having the potential bifurcation value x k in (4-7), we now prove in the following 
theorem that the steady state bifurcation occurs at ( u, v, w, x k ) for each k G N + , which 
establishes existence of nonconstant positive solutions to (4.1). 


Theorem 4.2. Assume that a 3 > An + f3 32 and o(w) < 0. Suppose that for positive 
integers k,j G N + , 


X S k A Xj, Vfc A j and x k A Xk , VA: G N + , 


(4.11) 
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where x k and x k are given by (3.4) and (3.5) respectively. Then for each k G N + , there 

exist a positive constants and a unique one-parameter curve T k (s) = {(u k (s, x),v k (s, x),w k (s, x),Xk(s)) : 

s G (—5, 5)} of spatially inhomogeneous solutions ( u, v,w,x) G X x X x X x M to (4.1) 

that bifurcate from (u, v, w) at x = Xk- Moreover, the solutions are smooth functions of s 

such that 


Xk(s ) = Xk + °( s ), s e 5), 


(4.12) 


and 


(■ u k (s,x),v k (s,x),w k (s,x )) = (u,v,w) + s(u k ,v k ,w k ) + 0(s 2 ),s G (-5,5), (4.13) 


where ( u k ,v k ,w k ) is given by (4.8) and 0(s' 2 ) G Z is in the closed complement of 
J\f(D(F(u, v, w, x)) defined by 


Z = j(w, v, w) G X x X x X J 


uu k + vv k + ww k dx — 0 >. 


(4.14) 


Proof. All the necessary conditions except the following have been verified in order to 
apply the Crandall-Rabinowitz local theory in [9] 


v, w, x))(u k , v k , w k )\ x=x s 0 n{DF{u, v, w, %)), 


(4.15) 


where 1Z is the range of the operator. We argue by contradiction and suppose that condition 
(4.15) fails, then there exists a nontrivial solution (u, v, w) that satisfies 


' d\u" — x k il(j>(w)w" — a\uu + f3\uw = (- : £) 2 u(j)(w ) cos x G (0, L), 

d 2 v" — £ V(f>{w)w" — a 2 vv + (3 2 vw = 0, x G (0, L), 

d 3 w" — /3 3 i wu — (3 32 wv — a 3 ww = 0, x G (0, L). 

^u'(x) = v'(x) = w'(x) = 0, x — 0, L. 


(4.16) 


Multiplying equations in (4.16) by cos^ and integrating them over (0, L) by parts, we 
obtain that 




«i u 


0 

~/3 3 iw 

' (kTr) 2 u<t>(w )' 

2 L 
0 

0 


0 Xk^Hw)^) 2 +/3iu\ (f 0 ucos^dx\ 

-d 2 (!f ) 2 - a 2 v £A/>(u ;)(^) 2 + fi 2 v I v cos ^ffdx 
-fi 32 w -d 3 (!f) 2 -a 3 w J [f°w cos ^dxj 

(4.17) 


The coefficient matrix is singular because of (4.7), then we reach a contradiction and this 
completes the proof of condition (4.15). Finally the statements in Theorem 4.2 follow 
from Theorem 1.7 of [9]. □ 
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4.2 Stability of bifurcating solutions near (u,v,w,Xk) 

Now we proceed to study the stability of the spatially inhomogeneous solution (u k (s, x),v k (s, x), 
w k (s, x)) established in Theorem 4.2. Here the stability or instability is that of the bifurca¬ 
tion solution regarded as an equilibrium of system (3.1). To this end, we want to determine 
the turning direction of the bifurcation branch T k (s) around each bifurcation point yf. It 
is easy to see that the operator T is C 4 -smooth if 0 is C' 5 -smooth, therefore according to 
Theorem 1.8 in [ 10], we can write the following expansions 

'u k (s, x) — U + sP k COS + S 2 Lpi(x) + S 3 (p 2 (x) + o(s 3 ), 

Vk(s, x) = V + sQ k cos + s 2 0i(x) + s 3 0 2 (x) + o(s 3 ), ^ 18 ^ 

w k (s,x) = U! + s cos ^ + s 2 7 i(x) + s 3 y 2 (x) + o(s 3 ), 

,Xk(s) = Xk + S ^1 + s2 fc 2 + °( s2 ) ; 


where (</?*, 70 e Z in (4.14) and 1C, are constants for i — 1,2. o(s 3 ) are taken with 

respect to the T'-topology and o(s 2 ) is a constant. Moreover we have the Taylor expansion 


<j)(Wk(s, x)) = (j>(w) + S(f)'{w ) cos 


k'KX 

~T~ 



+ o( S 2 ). 

(4.19) 


First of all, we claim that the bifurcation branch T*. is of pitch-fork type by showing 
/Ci = 0. Substituting (4.18) into (4.1) and collecting s 2 terms, we obtain the following 
system 


ditpi - xfiZ0(w)7i = {aiipi - f3iyi)u - /Ci?20(u))(^) 2 cos ^ + Rk, x G (0 ,L), 

d 2 ip‘l - £v(j)(w)'Y , { = (a 2 0i - /3 2 7 i )v + S k , x <E (0, L), 

^37l = (031-Pfe + 032<3fc + 03 ) COS 2 ^ + (031 V?1 + 03201 + «37l)« ) > X G (0, L), 

= 0l(x) = 7i(x) = 0, x = 0,L, 

(4.20) 


where 




-Xk{-^) 2 i U( l>'( w ) + p kH w )) cos — 


0 l)Pfc cos 2 


IT’ 


and 




- 0 (^) 2 (v 0 '(w) + Qk<f>(w)) cos + (a 2 Qfc 


0 2 )Q fc cos 2 


knx 

~T 


Multiplying the equations in (4.20) by cos^f and integrating them over (0, L ) by parts 
give us 



(4.21) 
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(d 2 (-^-) 2 + a 2 v)^j J i>! cos ~^dx+ t,v(j)(w)(^-) 2 - /3 2 v j J 71 cos -j~dx = 0 , 

(4.22) 


and 

n - f L toTX _ _ f L , k-KX , / , Jc 7 T 2 _\ f L kirx 

P31VJ (pi cos —— dx + P3 2 w ipi cos —— ax + ( d3( —) +a3Wj 71 cos 


/ 0 


L 


1 0 


L 


L 


0 


L 


Since (<pi, 0i, 71 ) G Z, we infer from (4.14) that 


r^ /c7r x 

P k / ip 1 cos -j—dx + Q k 


f L , /c7ra /‘ L 

/ -01 cos ——ax + / 71 cos ——ax = 0 , 

'0 L Jo L 


dx = 0. 
(4.23) 


(4.24) 


where P^, Qfc are given by (4.9) and (4.10). Combining (4.22)-(4.24) gives 

0 d 2 ( ! f) 2 + a 2 v -£v0(w )(^) 2 - /3 2 v\ f pi cos ^ffdx^ 

fulfil cos dx 
J 0 L 7l cos ^dxj 


P 31 W 

p k 


P32W 

Qk 


d3(f) 2 + a3W 


The determinant of the coefficient matrix At to the system above is 


det(Al) =Pk(d 2 (^-) 2 + a 2 v^j (d 3 (^r) 2 + a 3 v?) + Q k Pziw(^ - £v<j>(w)(^-) 


.kit 


kn s 


- /? 3 i w(d 2 (^-) 2 + a 2 v^j - P k /3 32 w(^ - £n0(w)(^) 
=P k {H 2 H3 + P 3 2wQ k H 2 ) - p 31 w(Q 2 k H 2 + H 2 ) 


kn . 


/?32 ~ P.i 

/?31 

(f3 32 wH 2 


— ( — “5 Qk — 7T {H2H3 + fi 3 2 wQ k H 2 ) — /33iw(Q 2 H 2 + P 2 

P31 


in 2 p — zj \ rP2. (H 2 H 3 P -zj\ '^32^2^3^ 

-h 03\WH 2 jQ k — ^ — — + (33iwH 2 J - — -Qfc < 0. 


2 / 332 ^ 2-^3 , 


v & 

Therefore we have 


<y7l COS 


knx 


, f kirx , f kirx , .. 

ax = cos ——dx = / 71 cos ——ax = 0, (4.25) 

./n -C/ ./n -C/ 


and this implies that /Ci = 0 in (4.21). Thus the bifurcation branch Tps) is of pitch-fork, 
i.e., being one sided. Now we present another main result of this paper which states that 
the stability of the bifurcating solutions depends on the sign of /C 2 . 


Theorem 4.3. Suppose that all the conditions in Theorem 4.2 are satisfied and let Y k (s) = 
{(u k (s, k),v k (s, k),w k (s, k),Xk(s))} be the bifurcation branch given by (4.12)-(4.13). 
Denote Xo = mm fceN +{x£, Xk } as in (3.3). Then it holds that: (i) If Xo = Xk 0 < 
min fceN + x k , then Y ko (s) around (u,v,w,x f 0 ) is asymptotically stable when /C 2 > 0 
and it is unstable when /C 2 < 0, while Y k (s) around (u,v,w,x k ) is always unstable for 
each k f ko; (ii) If Xo = X kl < min fceN + x k , then Y k (s) around (u,v,w,x k ) is always 
unstable for each k G N + . 
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The bifurcation curves I\.(s) in case (i) are illustrated in Figure 1 schematically. Our 
results suggest that (u,v,w,x k ) loses its stability to stable steady state bifurcating so¬ 
lution with wave mode number for which y k achieves its minimum over N + . When 
case (ii) occurs, we surmise that the stability of the homogeneous solution is lost to sta¬ 
ble Hopf bifurcating solutions. This is rigorously verified in Section 5. We would like to 
mention that, /C 2 can be evaluated in terms of system parameters and we give the detailed 
calculations in Appendix. 




Figure 1: Pitch-fork bifurcation diagrams when case (i) in Theorem 4.3 occurs. The stable 
bifurcation curve is plotted in solid line and the unstable bifurcation curve is plotted in 
dashed line. The branch r ko (s) around (u, v, w, x ko ) is stable if it turns to the right and 
is unstable if it turns to the left, while F/,.(s) around (u, v, tv, yf.) is always unstable if 

k k 0 . 


Proof of Theorem 4.3. Our proof follows the approaches in [5i , 6 ] based on slight mod¬ 
ifications in the arguments for Corollary 1.13 of [ 0], or Theorem 3.2 of [5' ], Theorem 
5.5, Theorem 5.6 of [8]. We shall only prove case (ii) and case (i) can be treated similarly. 

For each k G N + , we linearize (4.1) around (u k (s, x), v k (s, x), w k (s, x), y fc (s)) and 
obtain the following eigenvalue problem 

DJ r (u k (s,x),v k (s,x),w k (s,x),Xk(s))(u,v,w) = cr(s)(u,v,w), ( u,v,w ) G TxTxT, 

then (u k (s, x), v k (s, x),w k (s, x), Xk{$)) is asymptotically stable if and only if the real part 
of eigenvalue o(s) is negative. 

Sending s —> 0, we know from the proof of (4.15) that o = 0 is a simple eigenvalue 
of DIF(u , v, w , xi ) = v i w ) or equivalently 

'd\u" — x k uf(w)w" — ot\UU + f3iuw = cm, x G (0, L), 

d 2 v" — £,v(f)(w)w" — a 2 vv + (3 2 vw = ov, x G (0, L ), 

d^w" — /? 3 i vju — / 33 2 wv — QL 3 WW = aw, x G (0, L). 

^u'(x) = v'(x) = w'(x) =0, x — 0, L, 
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which has one-dimensional eigen-space Af(DT(u , v, w, xf)) = {(-Pfc, Qk, 1) cos 
Multiplying the system above by cos^f and integrating them over (0, L) by parts, we 
have that o = 0 is an eigenvalue of (3.2) with x = xf, which reads 



If xo = min fceN + Xk < Xk for a11 k e or Xo = min fceN + xl < Xk for k + k o, we 

have from the proof of Proposition 3.1 that this matrix always has an eigenvalue a with 
positive real part. From the standard eigenvalue perturbation theory in [ ], for s being 

small, there exists an eigenvalue cr(s) to the linearized problem above that has a positive 
real part and therefore (uk(s,x),Vk(s,x),Wk(s,x),Xk(s)) is unstable for s G (—5,5). □ 

According to Theorem 4.3, the only stable bifurcation branch must be Tf (s) if Xk 0 — 
min feeN +{xf, xjf}> therefore ( u,v,w ) loses its stability only to nonconstant steady state 
with wave mode COS This gives a wave mode selection mechanism for system (1.1) 
when x is around the bifurcation value. In general it is very difficult to determine whether 
Xo is achieved at xf or xjf • According to the discussions after Remark 3.1, if the interval 
length L is sufficiently small, xo = xf < niin / r. eN + Xk an d the only stable bifurcating 
solution has wave mode cos p which is spatially monotone. The wave mode section 
mechanism given in Theorem 4.3 is verified and illustrated in our numerical studies of 
(1.1) in Section 6. 

5 Time-periodic positive solutions 

In this section, we study the periodic orbits of (3.1) that bifurcate from (u,v,w) at x = 
Xk ■ We want to show that under proper assumptions on system parameters, the constant 
equilibrium (u, v. w) loses its stability through Hopf bifurcation as x comes across xo = 
min fceN +{xf, Xk}- To apply the bifurcation theory for (3.1) at point x = Xjf> we need to 
verify that the real part of eigenvalue crosses the imaginary axis at Xk- 

According to the discussions in Section 3, Hopf bifurcation occurs for (3.1) at (u, v, w) 
only if x = Xk an< ^ 7 /i (X- k ) > 0, when the stability matrix (3.2) has purely imaginary 
eigenvalues given by 



To determine when rft (xjf , k) > 0, we let xt 1 be the unique root of ?/i(x, k) — 0 which is 
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given explicitly in the following form 


X k 


M 


+ a 1 u)(d- 2 ( ! f) 2 + a 2 v) + (5i(^) 2 + aiu)(d 3 (^) 2 + a 3 w) 


P31 uwcpiw )^) 2 

(4(x) 2 + «2^)(^3(x) 2 + a sw) C@32V(j){w){^) 2 + @ 32 @ 2 V + @31@lU 


@31 uw<j)(w)(!f) 2 


@31 «*(«0® 2 


We first give the following fact which will be used in our coming analysis. 

Lemma 5.1. Let \jl l be given as above, then for each k 6 N + ,it holds that either xj r < 


Xk < Xk or Xk < Xk < Xk- Moreover rji(x%, k) > 0 > 77i(xf, k) ifXk < Xk < Xk 
and /7i(Xfe, k) > 0 > m(Xk> k ) l fXk < Xk < Xk- 


According to Lemma 5.1 and discussions in Section 3, the stability matrix (3.2) has 
a pair of purely imaginary eigenvalues if and only if x — Xk < Xh> therefore Hopf 
bifurcation may occur at (u, v, w, Xk) only when \k < Xk- We shall always assume this 
condition in the coming Hopf bifurcation analysis. 

5.1 Hopf bifurcation 

In this subsection, we prove the existence of Hopf bifurcation of (3.1) assuming that Xk < 
Xk- We recall the notation of Sobolev space X = {u G H 2 (0, L)\u'(0) = u'(L ) = 0} 
from Section 4. According to the proof of Theorem 2.1, we know that (3.1) is normally 
parabolic, therefore we can apply the Hopf bifurcation theory from [5] (or Theorem 1.11 
from [1 ], Theorem 6.1 from [36]). Our main result on the existence of nontrivial periodic 
orbits of (3.1) states as follows. 

Theorem 5 . 2 . Suppose that all parameters in (3.1) are positive, 03 > /j 3 i + T32 and 
<j)(w) < 0. Assume that x] f 7^ xf f or Y) f k and Xk < Xk’ then there exist a posi¬ 
tive constant 5 and a unique one-parameter family of nontrivial periodic orbits Okif = 



such that ( u k (s, x, t ), Xk(s )) is a nontrivial solution of (3.1) and Uk{s, x, t ) is periodic of 
time t with period 



(5.2) 


and {(V^, ±iTo)} cire eigen-pairs of matrix (3.2); moreover i)k(si) f dk(s 2 ) for all 
-Si f s 2 G (—<), 5) and all nontrivial periodic solutions around (u, v, w, Xk) must be on 
the orbit r 0k{s)- s e (—5,5). In other words, if (3.1) has a nontrivial periodic solution 
u i(x, t) with period T for some x £ around dk(s) and a small positive constant e such 
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that \x — xf(*)l < e,\T — ^r \ < e and max (Et+ xe ^ \ui(x, t ) — ( u , v,w) \ < e, then there 
exist constants s 0 G {—5,5) and 9 0 G [0, 2n) such that (T, y) = (T fc (s 0 ), X k { s o)) an d 
Ui(x,t) = u k {s 0 ,x,t + 9 0 ). 


Proof. We follow the approach in the proof of Theorem 5.2 in [59] or Theorem 3.4 in [36]. 
According to Proposition 3.1 and Remark 3.1, the stability matrix (3.2) with x — X k has 
a pair of purely imaginary eigenvalues crff 3 (fc) = ± y///i {x k )f moreover since x k 7^ xf 
for Vj f k, matrix (3.2) has no eigenvalue of the form mToi for m G N + \{±1}. 

Let (Ji{Xi k) and crf^fy, k) = A(y, k ) ± ir{x, k) be the unique eigenvalues of (3.2) in 
a neighbourhood of x = X k ■ Then erf 7 , A and r are real analytical functions of x satisfying 
A {x k , k) — 0 and r{x k , k) = r 0 > 0. In order to apply Hopf bifurcation theory, we need 
to prove the following transversality condition 


d\{x, k) 


dx 


X=Xk 


7^0. 


(5.3) 


Substituting the eigenvalues erf^y, k) and cr^f x, k) = A(y, k) ± ir{x, k) into the char¬ 
acteristic equation of the stability matrix (3.2) and equating the real and imaginary parts 
give 

[ ~V 2 (x,k) = 2\(x,k) + o?(x,k), 

{ Vi(x,h) = X 2 (x,k) + r 2 (x,k) + 2\{x,k)a?{x,k), (5.4) 

{ -Vo(x, k ) = (A 2 C\, k ) + r 2 (y, k))oi {x, k). 

Differentiating the equations above with respect to y, we obtain 


2A'(y,A;) + cri ix,k) = 0 


and 


2A(y, k) A'(y, k) + 2r(y, fc)r'(y, k) + 2A'(y, k)a x {x, k) + 2A(y, k)a[(x, k) 

=An uwfiw)^) 2 , 

(2A(y, k) A'(y, k) + 2r(x, k)r'{x, k))<n(x, k ) + (A 2 (x, k ) + ^ 2 (x, &)K(x, *0 
= - f3 3 iuwf{w){^) 2 {d 2 {^) 2 + a 2 v)- (5.5) 

Since A(yf, fc) = 0 and cri(xf, fc) = -r) 2 {x k , solving (5.5) with y = xf gives that 


_ /3 3 iw^(w)(^) 2 (d 2 (^) 2 + q 2 t> - r/ 2 (yf, fc)) 

r o + vZUk, k) 

fi3iWU(f){w) {kjf) 2 {{d\ + d 3 )(x) 2 + "1“ + a 3^) 

r o + ^(xf, A:) 


and 

A'(xf,fc) = -^i(xf) > 0. 

This verifies all the transversality conditions required in applying the Hopf bifurcation 
theory, then Theorem 5.2 follows from Theorem 1 in [5]. □ 
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Theorem 5.2 implies that system (3.1) admits time-periodic spatial patterns that bi¬ 
furcate from (u,v,w,Xk) if and on ly if xl! < xf • Furthermore, it gives the explicit 
expression of the time-periodic spatial patterns as 0 k mentioned above with the spatial 
profile of eigen-function cos pp 

As we have discussed in Section 3, it is very difficult to determine the necessary con¬ 
dition Xk < Xk i n tcrms °f system parameters, however if the interval L is sufficiently 
small, we always have x k > Xk f° r cac h k G N + , and therefore this indicates that Hopf 
bifurcation dose not occur for (3.1) when the interval length is sufficiently small. Indeed, 
in this case, we already know from the discussions after the proof of Theorem 4.3 that the 
stability of the homogeneous solution (u, v. w) is lost through the steady state bifurcation 
at the first bifurcation branch (■ u,v,w , xf )> which contains stable stationary solutions of 
(3.1) with eigenfunction cos(p). 


5.2 Stability of time-periodic bifurcating solutions 


We continue to explore the stability of the time-periodic bifurcating solutions on the bi¬ 
furcation curves d k (s) obtained in Theorem 5.2. The stability here we mean is the formal 
linearized stability of a periodic solution relative to perturbations from d k (s). Suppose 
that Xk 0 = min fceN + x k < xf, V/c G N + , and assume that all the conditions in Theorem 
5.2 are satisfied here, then our stability results show that d k (s), s G (—5, 5) is asymptoti¬ 
cally stable only if x — X ko - 

Denote u k (s,x,t) = (u k (s,x,t),v k (s,x,t),w k (s,x,t)) and let (u k (s,x,t),T k (s), 
Xk( s )) be the periodic solutions on the branch d k (s) obtained in Theorem 5.2. Then we 
can rewrite (3.1) into the following form 

= £(ufc,Xfc(s)), 

where 

/ {dxu' k - Xk{s)uk^{w k )w' k y + «i(l - u k )u k + /3i u k w k 
Ufc, Xfc(s)) = (d 2 v k - ivk^Wkjw'b)' + a 2 (l - v k )v k + /3 2 v k w k 
\ dyw" + a 3 (l - w k )w k - foiu k w k - /3 32 v k w k 



Differentiating the system against t, writing ii k — we have 


du k 

dt 




then we observe that 0 is a Floquet exponent and 1 is a Floquet multiplier for u^. 

Linearize the periodic solution around the bifurcation branch d k (s) by substituting the 
perturbed solution + we~ lt , where w is a sufficiently small T-periodic function and 
l — l(s) is a continuous function of s, then we have that 


dw(s,t) 


dt 


G u {»k, Xk(s))vr(s, t ) + /(s)w(s, t ), 


(5.7) 
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where Q u is the Frechet derivative with respect to u given by 

Sll(ufc, Xfc(s)) = Dg(u k ,v k ,w k ,Xk(s))(u; V, w ) 

( d]_u" - Xk( a )(. u <K™k)™' k + u k wc/>'(w k )w k + u k tl>(w k )w'y + (ai - 2a 1 u k + pxw k )u + /3iUfcuA 
d 2 v" - i(vt/>(w k )ui' k + v k wt/>'(w k )w' k + v k c/>(w k )w')' + (a 2 - 2a 2 v k + p 2 ™k) v + p 2 V k ™ ) ■ 
d 3 w" - p 31 u k w - p 32 v k w + (q 3 - 2a 3 w k - fel«t - P 32 v k )w ] 

The stability of the bifurcating solutions around Xk can determined by computing the 
eigenvalues of this reduced equation. When s = 0, (5.7) is associated with the eigenvalue 
problem 

Qo(k)w = /(0)w, (5.8) 

where 

M o 0 ~XkU(j){w) 

Go(k) = I 0 d 2 £ 2 - a 2 v ~^v<j)(w)j^ +/3 2 v 1, 

V -031W -P 32 W -d^-a^w J 

the spectrum of which is infinitely dimensional. Moreover Q 0 corresponds to the stability 
matrix (3.2). 

(-di( J f ) 2 ~ aiu 0 xf #(^)(^) 2 + Piu\ 

A'(xf) = 0 -d 2 (f) 2 -a 2 v £vcj)(w)( J f ) 2 + (3 2 v \ ,j G N + . 

\ -foiw -/3 3 2W 4(x ) 2 _ «3 W / 

(5.9) 

Suppose that min feeN +{xjf, xf} = Xk 0 * or some k 0 e N + . We first show that i? fc (s) 
around Xk ’ s unstable for any k 7 ^ k 0 . Denote the eigenvalues of Ak{Xk) by a i(',Xk 1 k), 
a 2 (Xkik) and a^(Xkik). According to the Proposition 3.1, there exists at least one 
eigenvalue with positive real part if x > To- Therefore for any positive integer k =A A’ 0 , 
we have that Q 0 (k ) must have an eigenvalue with positive real part hence Z( 0 ) < 0 if k 7 ^ 
k 0 . By the standard perturbation theory for an eigenvalue of finite multiplicity [ , ], 

l(s) < 0 for s being small if k k 0 , therefore all the bifurcation branches t)/,(.S') around 
(u. v. w) are unstable if k ^ /c 0 . The result indicates that if a periodic bifurcation solution 
is stable, it must be on the $ fco ( s ) branch where Xk 0 < min feeN + xf> be., it is on the 
left-most branch, while the branches on its right hand side are always unstable. 

Now we proceed to discuss the stability of branch dk 0 (s) around ( u,v,w,Xk 0 )• Ac¬ 
cording to Lemma 2.10 in [11] (or [2^ , 25]), the eigenvalue l(k) is a continuous real 
function of s near the origin. For x being around the eigenvalue of (5.9) are <7i(x, k) 
and cr 2 ) 3 (x, k) = A(x, k) ± ir(x, k). According to Theorem 2.13 in [ ], l(s) and sx' ko (s) 

have the same zeros in small neighbourhood of s = 0 where l{k) and —A (Xk 0 ) s x'k 0 s ( s ) 
have the same sign — A(x^)sXfc 0 s(s) 7 ^ 0 , l(s) 7 ^ 0 , and 

\K S ) + A((xf 0 )sXfe 0 s(s)| < |sxi 0 (s)|o(l), as s -> 0. 


According to Theorem 8.2.3 in [ ], if l(s) > 0, the periodic bifurcation solutions 

are orbitally asymptotically stable and, if l(s) < 0 , the periodic bifurcation solutions 
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are orbitally unstable. We have proved that A (Xk 0 ) < 0 and l{s) and (s) has the same 
sign. Therefore, assuming that ;\ " o (0) p 0, if the branching solutions appear supercritical, 
they are stable and if they appear subcritical, they are unstable. Therefore, one need to 
compute x'k 0 and/or x" n similarly as in Section 4. The calculations are straightforward but 
complicated and we skip them here for simplicity. 


6 Numerical simulations 

This section is devoted to the numerical studies of system (3.1). We are motivated to in¬ 
vestigate the effects of prey-taxis on the formation of nontrivial patterns to this system. 
In particular, we show that the one-dimensional system admits both stationary and time- 
periodic solutions emerging from bifurcations. Moreover, we shall see that, when the 
prey-taxis rate x is taken to be greatly larger than the critical bifurcation value xo, (3-1) 
can develop various interesting patterns with striking structures such as spikes, propaga¬ 
tion, coarsening, etc. 


6.1 Stationary patterns 

In Table 1, we list the values of xf ar| d Xk given in (3.4) and (3.5) respectively. It is 
shown that their minimum value over N + is achieved at xt ~ 6.05, therefore {u. v. w) = 
(1.71,1.18, 0.71) loses its stability through steady state bifurcation with wave mode cos ^p 
In Figure 2, we plot numerical solutions of (3.1) subject to initial data (u 0 ,v 0 ,w 0 ) = 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Xk 

66.98 

18.98 

10.30 

7.49 

6.41 

6.05 

6.09 

6.36 

6.80 

Xk 

1204.2 

3550.8^ 

504.2( 

) 575. 8 ( 

) 705.5( 

OO 

r- 

OO 

o 

OO 

SO 

31336.9 

31619.1 


Table 1: Values of xf in (3.4) and Xk i n (3-5) for L = 7. The system parameters are chosen 
to be d\ = d 3 = 0.1, d 2 = 2, ol\ — f3\ — /3 2 — 0.5, a 2 = 2, a 3 = 1 and fd 3 \ = j3 32 = 0.1, 
£ = 0.5. The sensitivity function is oiw) = w( 0.1 — w) which models the group defense 
of the preys when population density surpasses 0.1. We see that min fceN +{xf, Xk} = Xe- 
Therefore (u, v. w) loses its stability to the stable wave mode COS ^p. This is numerically 
verified in Figure 2. 

(u, v , w) + (0.01, 0.01,0.01) cos 7 rx, which are small perturbations from the homogeneous 
equilibrium. We see that the initial data have spatial profiles in the form of cos irx, but the 
spatial-temporal patterns develop according to the stable wave mode cos pT. 

Numerical simulations in Figure 3 are devoted to verifying that (u. v, w) loses its sta¬ 
bility to steady state bifurcations when xo is achieved at xf 0 for this set of system param¬ 
eters. These simulations support the results on the stability of the bifurcating solutions in 
Theorem 4.3. 
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Figure 2: Formation of stationary patterns of (3.1) over = (0,7). System parame¬ 
ters in all the graphes are taken to be the same as in Table 1 except that y = 8, which 
is larger than yf ~ 6.05 given in Table 1. Initial data are (-u 0 , r>o, Wo) = (u,v,w) + 
(0.01, 0.01, 0.01) cos7rx, while the stable pattern has wave mode COS ^p. These graphes 
support our stability analysis of the bifurcating solutions. 


Interval length L 

1 

2 

3 

4 

5 

6 

7 

8 

k 0 

1 

2 

3 

4 

5 

5 

6 

7 

XO = Xkn 

6.09 

6.09 

6.09 

6.09 

6.09 

6.08 

6.05 

6.04 

Interval length L 

9 

10 

11 

12 

13 

14 

15 

16 

k 0 

8 

9 

10 

11 

12 

13 

14 

15 

XO = Xkn 

6.04 

6.03 

6.04 

6.04 

6.04 

6.04 

6.04 

6.04 


Table 2: Stable wave mode numbers and the corresponding bifurcation values Xo for dif¬ 
ferent interval lengthes. System parameters are chosen to be the same as those in Table 1. 
We see that the threshold value y 0 is always achieved at the steady state bifurcation point 
xf 0 - This table also indicates that larger intervals support higher wave modes. 

6.2 Time periodic patterns 

Our next set of numerical results are provided to demonstrate that (3.1) admits time- 
periodic patterns through Hopf bifurcations. To this end, we set system parameters to be 
d\ = g?3 = 1, g?2 = 0.1, ot\ = 0.02, o?2 = 0.04, ol 3 = 8 and j 3 \ = 0.05, /?2 = /?3i = @32 = 
0.5, while the sensitivity function is chosen <fi(w) = 0.05tu(0.2 — w). We shall show that 
the equilibrium (u,v,w) = (2.13,6.65,0.45) loses its stability to time-periodic orbits. 
Table 3 lists the values of xf- an d Xk when the interval length is L = 7. We see that the 
threshold value yo is achieved at X3 ~ 92.57. In Figure 4, we plot the numerical solutions 
of (3.1) subject to initial data ( u 0 ,v 0 ,w 0 ) = (u,v,vj) + (0.01,0.01, 0.01) cos ttx, small 
perturbations from the homogeneous equilibrium. The initial data have spatial profiles 
in the form of cos ttx, but the spatial-temporal patterns develop according to the stable 
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u(x,t) 
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Figure 3: Formation of stationary patterns of (3.1) over intervals with lengthes L — 9, 
11, 13 and 15. System parameters here are taken to be the same as those in Table 1 ex¬ 
cept that x = 8, which is slightly larger than Xk 0 ~ 6.04 given in Table 2. Initial data 
are ( u 0 ,v 0 ,w 0 ) = (u,v,w) + (0.01, 0.01, 0.01) cos irx. These graphes support our sta¬ 
bility analysis of the bifurcating solutions and indicate that large intervals support more 
aggregates than small intervals. 


k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Xk 

106.4 

98.63 

107.32 

122.53 

143.15 

169.05 

200.24 

236.80 

278.76 

Xk 

186.37 

96.73 

92.57 

105.46 

127.03 

155.24 

189.40 

229.24 

274.64 


Table 3: Values of Xk ' n (3-4) and Xk ' n (3-5) for L = 7. System parameters are d\ = 
d 3 = 1, d 2 = 0.01, = 0.02, « 2 = 0.04, a 3 = 8 and /3i = 0.05, /3 2 = f3 3 1 = /3 32 = 0.5, 

while the sensitivity function is 4>(w) = u>(0.2 — w ). We see that min feeN +{xf, Xk} — 
X 3 ■ Therefore (u. v, w) loses its stability to the time-periodic solutions with wave mode 
COS ^p. This is numerically verified in Figure 4. 

time-periodic patterns with wave mode cos ^p. 


Interval length L 

2 

3 
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5 

6 

7 

8 

9 

k 0 

1 

2 

3 

4 

5 

5 

6 

7 

XO = Xkn 

97.68 

92.13 

97.68 

91.49 

92.13 

92.57 

91.15 

92.13 

Interval length L 

10 

11 

12 

13 

14 

15 

16 

17 

k 0 

8 

9 

10 

11 

12 

13 

14 

15 

Xo = Xkn 

91.5 

91.2 

91.13 

91.21 

91.30 

91.49 

91.15 

91.40 


Table 4: Stable wave mode numbers and the corresponding bifurcation values xo for dif¬ 
ferent interval lengthes, where system parameters are chosen to be the same as those in 
Table 3. We see that the threshold value Xo is always achieved at the Hopf bifurcation 
Point Xk 0 - 

Numerical simulations in Figure 5 are devoted to verifying that the ( u,v,w ) loses its 
stability to Hopf bifurcations when xo is achieved at Xk 0 ’ where system parameters are 
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U(x,t) 



V(x,t) 



V(x,t) 



W(x,t) 




Figure 4: Formation of time-periodic spatial patterns of (3.1) over fl = (0, 7). System 
parameters in all the graphes are taken to be the same as in Table 3 except that x = 120, 
which is slightly larger than Xk 0 ~ 92.57 given in Table 3. Initial data are (w 0 , n 0 , w o) = 
(u, v , w) + (0.01, 0.01, 0.01) cos 7tx, however the stable oscillating patterns have spatial 
profile COS which emerge periodically. These plots support our stability analysis in 
Section 5. 


taken to be the same as those in Table 3. In particular, we select the interval lengths to be 
L = 9, 11, 13 and 15 respectively. These simulations support our results on the stability 
of the Hopf bifurcating solutions obtained in Theorem 5.2. 
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Figure 5: Formation of time-periodic spatial patterns of (3.1) over intervals with lengthes 
L = 9,11,13 and 15 respectively. System parameters in all the graphes are taken to be the 
same as in Table 3 except that x = 120, which is slightly larger than Xk Q given in Table 4. 
Initial data are (w 0 , Vo, w 0 ) = (u, v, w) + (0.01, 0.01, 0.01) cos7rx. These graphes support 
our stability analysis of the bifurcating solutions. 
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6.3 Other interesting patterns 

In Figure 6, we plot the formation of stable boundary spikes of (3.1) through traveling 
wave over Cl = (0, 7). System parameters are taken to be d\ = 5, d 2 = 0.5, d 3 = 1, 
«i = a 2 = 0.05, a 3 = 4, fdi = 0.3, /3 2 = 0.5, /3 31 = 1 ,/3 32 = §, £ = 0.1 and <p(w) = 
w(0.1 — w). Prey-taxis rate x — 3000 is greatly larger than the critical bifurcation value 
Xo = 956.79. We observe that the boundary spike is developed through traveling wave 
solution. However, rigorous analysis of qualitative properties of the propagating solutions 
is out of the scope of our paper. Finally, we present numerical simulations in Figure 7 to 


u(x,t) 



Space x 



Space x 


v(x,t) 



Space x 
v(x,t) 



Space x 


w(x,t) 



Space x 
w(x,t) 
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Figure 6 : Formation and development of boundary spike through wave propagation. Initial 
data are (u 0 ,v 0 ,w 0 ) = (pL,v,w) + (0.01, 0.01,0.01) cos nx. Prey-taxis rate x = 3000 
which is far away from the critical bifurcation value Xo = 956.79. 


show that when the prey-taxis is much larger than xo, (3.1) admits some other interesting 
and striking dynamics such as merging and emerging of spikes, irregular spatial-temporal 
oscillations etc. For example, Subplot (i) of Figure 7 shows that there occurs a coarsening 
process in (3.1) in which interior spikes of u(x,t ) shift to the boundary or the center 
to merge into another stable spike. We also observe the spontaneous emergence of stable 
interior spikes at time t « 100. All the parameters and the initial data in (3.1) are chosen to 
be the same as in Figure 2, except that x = 124, which is much far away from xo ~ 6.05. 
In subplot (ii), when the system parameters and the initial data are chosen to be the same 
as in Figure 4, except that x = 160, (3.1) initially develops time-periodic spatial patterns 
which are metastable. Then the oscillating patterns develop into stable stationary spikes. 
Time-periodic patterns and spontaneous initiation of interior multiple spikes are observed 
in subplots (iii) and (iv). 
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u(x,t) u(x,t) u(x,t) u(x,t) 



Space x Space x Space x Space x 


Figure 7: Pattern formations in (3.1) due to the effect of large prey-taxis rate Various 
interesting and complex spatial-temporal dynamics are observed in this system. 

7 Conclusions and discussions 

Our paper investigates population dynamics of a two-predator and one-prey model with 
prey-taxis, given by a 3 x 3 reaction-advection-diffusion system. It is proved that the sys¬ 
tem admits positive classical solution which is global and uniformly bounded in time over 
ID or 2D bounded domains. The same results are obtained for its parabolic-parabolic- 
elliptic counterpart for domains of arbitrary space dimension. 

Stability of the unique positive equilibrium is studied when the domain is a finite in¬ 
terval. It shows that both prey-taxis x an d sensitivity function 0 determines the linearized 
stability of this equilibrium. It is known (see [ ] e.g.) that, in contrast to chemotaxis 

[18, 19] or advection for competition system [5 ], prey-taxis stabilizes constant equilib¬ 
rium for one-predator and one-prey system. However, our result reveals that this is true 
only when there is no group defense in the preys, i.e., a huge amount of preys can ag¬ 
gregate and keep their predators away from the habitat. If the predators retreat from the 
habitat, which can be modeled by choosing 4 >(w) < 0, prey-taxis destabilizes the con¬ 
stant equilibrium, which becomes unstable as x surpasses xo — min fceN {xf, xf} given 
by (3.3). Therefore group defense is an important mechanism in the formation of nontriv¬ 
ial patterns in (1.1). 

We have obtained both stationary and time-periodic spatial patterns to the system over 
ID bounded interval through steady state bifurcation at x = Xk an d Hopf bifurcation 
a t x = Xk respectively. Stabilities of these bifurcating solutions are also investigated 
rigorously. It is proved that steady state bifurcation occurs at (u, v, w, xf) f° r each k G 
N + , however, only the /;: 0 -branch that turns to the right is stable if Xo = Xk 0 < Hiin feeW Xk ■ 
In other words, if the steady state bifurcation curve is stable, it must be on the left most 
branch on the %-axis. On the other hand, Hopf bifurcation (u , v, w, xl!) only if xl! < 
Xk , while only the left most branch can be stable. Moreover, our analysis indicates that 
small intervals only supports steady state bifurcations while large intervals may lead to 
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Hopf bifurcation when the system parameters are chosen properly. Extensive numerical 
simulations are performed to illustrate and support our theoretical findings. Apparently, 
the formation of these nontrivial patterns is due to the effect of large prey-taxis and prey 
group defense effect. 

Global existence and bounded are obtained for (1.1) over 2D and it is interesting to ask 
the same question for the system over higher dimensions. Logistic decays in the kinetics 
help to prevent finite or infinite time blow-ups, however, whether or not they are sufficient 
over higher dimensions, in particular when the prey-taxis rate is large, is unknown in the 
literature. 

Our bifurcation analysis is based on the local versions in [5, 9] etc. From the viewpoint 
of mathematical analysis, it is interesting to investigate the behavior or shape of these lo¬ 
cal branches, in particular in the study of positive steady states when large prey-taxis 
may lead to striking structures such as spikes and layers, etc. For example, according to 
the global theory of Rabinowitz [ ] and its developed version in [5( ], global continuum 

of Tfc(s) either intersects with the ;\-axis at another bifurcating point, or extends to in¬ 
finity, or intersects with a singular point. Populations growth terms in (3.1) inhibits the 
application of topology argument developed in [8, 61 ] etc. 

When the prey-taxis rate is around bifurcation values, our results provide almost a 
complete understanding of the spatial-temporal dynamics of (1.1) over ID. Further re¬ 
search is needed on its pattern formations when x is away from xo and in particular when 
it is sufficiently large. For example, rigorous analysis of the profile of the spikes obtained 
in numerical simulations can be an interesting problem to probe in the future. There are 
also some interesting problems such as the investigation of chaotic dynamics in ( 1 . 1 ) or 
bifurcation analysis of (1.1) over higher dimensions. It is also meaningful to ask about the 
biologically realistic traveling wave solutions to (3.1), compared to those for the system 
without prey-taxis obtained in [35]. 


8 Appendix 


This section is devoted to evaluating /C 2 in terms of system parameters in (4.1). We know 
from Theorem 5.2 and Theorem 4.3 that /Ci = 0 in (4.18) and the steady state bifurca¬ 
tion branch T fc (s) is pitch-fork; moreover /C 2 determines the turning direction hence the 
stability of the steady state bifurcation Tfc(s) around xt■ F° r the generality, we obtain the 
general expression of /C 2 for each branch r fc (s), k G N + . 
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We collect and equate 0(s :i ) terms in (4.1) through (4.18) to have that 

r di </?2 = XkU(j)(w )72 - /C 2 U(j)(w)( ! f) 2 cos + (ai<^ 2 - /3 i7 2 )w 
+ ^( 2 «iP fc - /3i)v?i - /3iPfc7i) cos + XaA, 
d' 2 ^P '2 = £v(j)(w )72 + ((2a 2 Q fc - /3 2 )Pl - /3 2 <3fc7i) cos 
< +(a 2 ^ 2 - /5 2 7 2 )n + £C 2 , ( 8 -l) 

<hH = (/^l + ^32'01 + (Ail Ac + fimQk + ^)7l^ cos 

+ (/^ 31 ( 4’2 + fd 32 ^ 2 C^ 3 )w, 

^(p' 2 (x) = ij)' 2 {x) = i 2 (x) =0, X = 0, L, 


where 


„ - ± ,,-s,kn . 2 k 7 TX 

C 1 =-U(j) (w)(—) 7i cos 

i kn 2 kirx 
- <p{w)(—) ipi cos 


'(w) + Pfc 0 (w)) (^)7i sin 


knx 


l' t *~~~ l + 

knx 


(u<j)'(w) 
kn . 0 /c7rx 


+ Pk 4 >(w) ^7i cos + (j)"(w)w + 2 P k (j)'(w)J (-^) 2 sin 2 cos 
- ^ 0 "(w)w + Pfc 0 '(tD)j(^) 2 COS 


/c7rx 


knx 

~T~ ’ 


and 


C 2 = - vct)'(w){^-) 2 7i cos - ( 2 x 0 '(w) + Q k 4 >(w)J (^)7^ sin 


fc7rx 


....k'K. 2( knx .f-s/kn.,. . knx 
n w ){-^) Wicos— - <p[w){— )V' 1 sm— + 


tx 


fc7rx fc7rx 


+ Qk4>(w) ^ 7 " cos + [4>”{w)w + 2 Qk4>'(w)^j (“r ) 2 sin2 cos 
^0"(tD)w; + Q k <p'{w) ) (^) 2 cos 


/c7rx 


We multiply the first equation in (8.1) by cos^p and integrate it over (0, L), which im¬ 


plies that 

K, 2 u(p(w)(kn ) 2 
2 L 


= [di(-^) 2 + ai« 


<£2 cos ^^dx - [xkU^{w ){^-) 2 + /?!«) 


fc7T N 


ao 

rL 2 knx 
71 cos —-— 


7 2 cos ^dx + ^aiPfc - ^/A + ^Xfc0M(^r) 2 ) 

<Pi cos <2 ^Ldx - ^ [x k u(t)'{w){^-) 2 + /A Ac + XfcA0(w)) 

/ 11 A/7T \ C ^ 

dx + (^aiP fc - -/A - -Xfc0(w)(—) 2 J J <pidx 
71 dx. 


kn ' rL 


'(«>)(—) 2 + ftp* 


( 8 . 2 ) 
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On the other hand, we test the second and the third equations in (8.1) by same method to 
obtain 


0 = {d 2 {^-) 2 + a 2 vj I 0 2 cos ‘^j^dx - ^m0(w)(^) 2 + fd 2 v^j 


kn , 


knx 

-17 ' 

knx , / 1 1 (/ _wA: 7 r. 2 \ 

72 cos— cte + (^a 2 Q fc - 2^2 + ) J 

0i cos ^-^dx - ^ i^vcp'[w )(^-) 2 + / 3 2 Qfc + iQkfyiw)} 


(8.3) 


7i cos + ^a 2 Q fc - ^/ 3 2 - ^0(w)(^) 2 j / 0ioh 


kir , 


^(£u 0 '(w )(^) 2 +/3 2 Q fe ) / 7 idx, 


and 


_ „ _ f L knx , _ f L kirx , , Jot 2 f L 

0 =p 3 iW / p 2 cos——dx +/J 32 w / 0 2 cos —— cte + d 3 ( —) / 


+ 


/?31 


cos 


L 
2knx 


72 cos 


/otx 




L * + T 


"01 COS 


2kTTX 


dx 


'0 


+ 7(^31 Pk + PmQk + -77^, 

2 A 3 

/?32 


rL 2fc7ra 0 3 i /‘ l 

71 cos —-—dx H—— / (/9idx 


'0 


'0 


2 a.-I 


H— T - / il>idx + -(/3 3 i Pk + fd 32 Qk + 7~^) / 7i ( ^ x 
2 In 2 K-i In 


Moreover, from (<^ 2 , 02,72) £ -2 defined in (4.14), it follows that 


Pk 


kirx 


ki\x 


p 2 cos ——dx + Qk 02 cos ——dx + 


72 COS 


Aoto; 


dx = 0 . 


(8.4) 


(8.5) 


We combine (8.3)—(8.5) in the following system 

0 d 2 (!f) 2 + a 2 v -^fi 0 (w )(^) 2 - 0 ifi\ f <p 2 cos h f 1 dx\ 

foiw 032^ d 3 (^) 2 + a 3 w rL 

Pk Qk 1 


/7c 

Jq 02 COS ^jfdx 
Jo L 72 cos ^ dx) 


where 


Mi = - (a 2 <5fe - ^/?2 + ^ 0 (^)(^) 2 ) / 0 


2k'xx 


1 cos 


+ ^(£W>'M (^) 2 + 02<5fc + / 71 COS ‘^^dx 


dx 


2knx 



( 8 . 6 ) 
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( ' oi 2 Qk 


ft - ft(®xft 2 


M; )(“^“) 2 + fiiQk 


ft dx 


(8.7) 


7i dx , 


and 
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2kirx ft 2 
<f i cos —-—dx -— 


" L 2knx 

ft cos —-—dx 
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'0 
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Solving (8.6) by Cramer’s rule, we obtain that 
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' kn ^2 


ftft 


ft(ft ) 2 + a 3 w 


( 8 . 12 ) 


Due to the Neumann boundary conditions and ft = 0, integrating (2.23) by parts 
yields 

{ «iw / 0 L ^idx - ftw f 0 L jidx = -§P k (a x P h - ft), 

« 2 h Jq ft da; - (3 2 v / 0 L 7 i dx = ~^Q k {oi 2 Q k ~ ft), 

ftiw / 0 L <pda; + ft 2 w f 0 L ft da; + a 3 w / 0 L 71 da; = -f (fti P k + /ftA + « 3 ), 
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or equivalently 
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Solving (8.13) leads us to 
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Multiplying (2.23) by cos and integrating by parts yields 
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We have the solutions of (8.17) that 
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where the notations are 
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and 

/ 4 d 1 (!f) 2 + ai u 0 -XkU^wy'ff - Piu\ 

C 0 = 0 4d 2 (^) 2 + « 2 n -^{w)^) 2 - fd 2 v . (8.21) 

\ foivo fd 32 w 4d 3 (f) 2 + a 3 w J 

Note that JC 2 terms in (8.2) consists of integrals (8.9)-(8.10), (8.14) and (8.18)-(8.19). 
Therefore, given all the system parameters, we will be able to evaluate /C 2 and determine 
the stability of T fco (s) thanks to Theorem 4.3. 


References 

[1] P. Abrams and H. Matsuda, Effects of adaptive predatory and anti-predator be¬ 
haviour in a two-prey-one-predator system. Evolutionary Ecology, 7 (1993), 312— 
326. 

[2] B. E. Ainseba, M. Bendahmane and A. Noussair, A reaction-diffusion system mod¬ 
eling predator-prey with prey-taxis. Nonlinear Anal. Real World Appl., 9 (2008), 
2086-2105. 

[3] N. Alikakos, L p bounds of solutions of reaction-diffusion equations. Comm. Partial 
Differential Equations, 4 (1979), 827-868. 





Patterns of systems with prey-taxis 


43 


[4] H. Amann, Dynamic theory of quasilinear parabolic equations. II. Reaction- 
diffusion systems, Differential Integral Equations, 3 (1990), 13-75. 

[5] H. Amann, Hopf bifurcation in quasilinear reaction-diffusion systems, Delay Dif¬ 
ferential Equations and Dynamical Systems, Lecture Notes in Mathematics, 1475 
(1991), 53-63. 

[6] H. Amann, Nonhomogeneous linear and quasilinear elliptic and parabolic bound¬ 
ary value problems. Function Spaces, differential operators and nonlinear Analysis, 
Teubner, Stuttgart, Leipzig, 133 (1993), 9-126. 

[7] A. Chakraborty, M. Singh, D. Lucy and P. Ridland, Predator-prey model with 
prey-taxis and diffusion. Math. Comput. Modelling, 46 (2007), 482-498. 

[8] A. Chertock, A. Kurganov, X. Wang and Y. Wu, On a chemotaxis model with 
saturated chemotactic flux, Kinet. Relat. Models, 5 (2012), 51-95. 

[9] M. G. Crandall and P. H. Rabinowitz, Bifurcation from simple eigenvalues, J. 
Functional Analysis, 8 (1971), 321-340. 

[10] M. G. Crandall and P. H. Rabinowitz, Bifurcation, perturbation of simple eigenval¬ 
ues, and linearized stability. Arch. Rational Mech. Anal., 52 (1973), 161-180. 

[11] M. G. Crandall and P. H. Rabinowitz, The Hopf bifurcation theorem in infinite 
dimensions. Arch. Rational Mech. Anal., 67 (1977), 53-72. 

[12] T. Czaran, Spatiotemporal Models of Population and Community Dynamics, Chap¬ 
man and Hall, London, 1998. 

[13] [10.1111/j.1469-1809.1937.tb02153.x] R. A. Fisher, The wave of advance of ad¬ 
vantageous genes. Annals of Eugenics, 7 (1937), 355-369. 

[14] D. Griinbaum, Advection-diffusion equations for generalized tactic searching be¬ 
haviours, J. Math. Biol., 38 (1999), 169-194. 

[15] [10.1016/0022-5193(71)90189-5] W. D. Hamilton, Geometry for the selfish herd, 
J. Theoret. Biol., 31 (1971), 295-311. 

[16] X. He and S. Zheng, Globed boundedness of solutions in a reaction-diffusion system 
of predator-prey model with prey-taxis, Appl. Math. Lett., 49 (2015), 73-77. 

[17] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Springer-Verlag- 
Berlin-New York, 1981. 



44 


Ke Wang, Qi Wang and Feng Yu 


[18] T. Hillen and K. J. Painter, A user’s guidence to PDE models for chemotaxis, J. 
Math. Biol., 58 (2009), 183-217. 

[19] D. Horstmann, 1970 until present: The Keller-Segel model in chemotaxis and its 
consequences. I, Jahresber DMV, 105 (2003), 103-165. 

[20] D. Horstmann, Generalizing the Keller-Segel model: Lyapunov functionals, steady 
state analysis, and blow-up results for multi-species chemotaxis models in the pres¬ 
ence of attraction and repulsion between competitive interacting species, J. Nonlin¬ 
ear Sci., 21 (2011), 231-270. 

[21] D. Horstmann and M. Winkler, Boundedness vs. blow-up in a chemotaxis system, 
J. Differential Equations, 215 (2005), 52-107. 

[22] L. Hsiao and P. de Mottoni, Persistence in reading-diffusing systems: Interaction 
of two predators and one prey. Nonlinear Anal., 11 (1987), 877-891. 

[23] L. Jin, Q. Wang and Z. Zhang, Qualitative Studies of Advective Com¬ 
petition System with Beddington-DeAngelis Functional Response, preprint, 
http: //arxi v. org/ab s/1412.3371 

[24] D. D. Joseph and D. Nield, Stability of bifurcating time-periodic and steady solu¬ 
tions of arbitrary amplitude. Arch. Rational Mech. Anal., 58 (1975), 369-380. 

[25] D. D. Joseph and D. H. Sattinger, Bifurcating time periodic solutions and their 
stability. Arch. Rational Mech. Anal., 45, (1972), 75-109. 

[26] [10.1086/284707] P. Kareiva and G. Odell, Swarms of predators exhibit ’’preytaxis” 
if individual predators use area-restricted search. The American Naturalist, 130 
(1987), 233-270. 

[27] T. Kato, Functional Analysis, Springer Classics in Mathematics, 1995. 

[28] [10.1016/0022-5193(70)90092-5] E. F. Keller and F. A. Segel, Initiation of slime 
mold aggregation viewed as an instability, J. Theoret. Biol., 26 (1970), 399-415. 

[29] K. Kuto, Stability of steady-state solutions to a prey-predator system with cross- 
diffusion, J. Differential Equations, 197 (2004), 293-314. 

[30] K. Kuto and Y. Yamada, Multiple coexistence states for a prey-predator system with 
cross-diffusion, J. Differential Equations, 197 (2004), 315-348. 

[31] O. A. Fadyzenskaja, V. A. Solonnikov and N. N. Ural’ceva, Finear and Quasi-Finear 
Equations of Parabolic Type, American Mathematical Society, 1968, 648 pages. 



Patterns of systems with prey-taxis 


45 


[32] J. M. Lee, T. Hilllen and M. A. Lewis, Continuous traveling waves for prey-taxis, 
Bull. Math. Biol., 70 (2008), 65^676. 

[33] J. M. Lee, T. Hilllen and M. A. Lewis, Pattern formation in prey-taxis systems, J. 
Biol. Dyn., 3 (2009), 551-573. 

[34] C. Li, X. Wang and Y. Shao, Steady states of a predator-prey model with prey-taxis, 
Nonlinear Anal., 97 (2014), 155-168. 

[35] J.-J. Lin, W. Wang, C. Zhao and T.-H. Yang Global dynamics and traveling wave 
solutions of two predators-one prey models. Discrete Contin. Dyn. Syst-Series B, 
20 (2015), 1135-1154. 

[36] P. Liu, J. Shi and Z.-A. Wang, Pattern formation of the attraction-repulsion Keller- 
Segel system. Discrete Contin. Dyn. Syst-Series B, 18 (2013), 2597-2625. 

[37] [10.1016/S0040-5809(03)00105-9] I. Loladze, Y. Kuang, J.-J. Elser and W.-F. Fa¬ 
gan, Competition and stoichiometry: Coexistence of two predators on one prey, 
Theoret. Pop. Biol., 65 (2004), 1-15. 

[38] [10.4319/lo.2013.58.5.1621] Z. Maciej Gliwicz, P. Maszczyk. J. Jabloriski and D. 
Wrzosek, Patch exploitation by planktivorous fish and the concept of aggregation as 
an antipredation defense in zooplankton, Fimnology and Oceanography, 58 (2013), 
1621-1639. 

[39] M. Mimura and K. Kawasaki, Spatial segregation in competitive interaction- 
diffusion equations, J. Math. Biol., 9 (1980), 49-64. 

[40] J. D. Murray, Mathematical Biology, Springer, New York, 1993. 

[41] T. Nagai, T. Senba and K. Yoshida, Application of the Trudinger-Moser inequality 
to a parabolic system of chemotaxis, Funkcial. Ekvac., 40 (1997), 411 —433. 

[42] W. Nagata and S.-M. Merchant, Wave train selection behind invasion fronts in 
reaction-diffusion predator-prey models, Phys. D, 239 (2010), 1670-1680. 

[43] K. Nakashima and Y. Yamada, Positive steady states for prey-predator models with 
cross-diffusion, Adv. Differential Equations, 1 (1996), 1099-1122. 

[44] A. Okubo and S. A. Levin, Diffusion and Ecological Problems, Modern Perspec¬ 
tives, 2nd Edition, Springer-Verlag, New York, 2001. 

[45] P. Pang and M. Wang, Strategy and stationary pattern in a three-species predator- 
prey model, J. Differential Equations, 200 (2004), 245-273. 



46 


Ke Wang, Qi Wang and Feng Yu 


[46] C. S. Patlak, Random walk with persistence and external bias. Bull. Math. Biophys., 
15(1953), 311-338. 

[47] P. Rabinowitz, Some global results for nonlinear eigenvalue problems, J. Functional 
Analysis, 7 (1971), 487-513. 

[48] K. Ryu and I. Ahn, Positive steady-states for two interacting species models with 
linear self-cross diffusions. Discrete Contin. Dyn. Syst., 9 (2003), 1049-1061. 

[49] [10.1086/375297] N. Sapoukhina, Y. Tyutyunov and R. Arditi, The role of prey 
taxis in biological control: A spatial theoretical model, The American Naturalist, 
162 (2003), 61-76. 

[50] J. Shi and X. Wang, On global bifurcation for quasilinear elliptic systems on 
bounded domains, J. Differential Equations, 246 (2009), 2788-2812. 

[51] N. Shigesada, K. Kawasaki and E. Teramoto, Spatial segregation of interacting 
species, J. Theoret. Biol., 79 (1979), 83-99. 

[52] G. Simonett, Center manifolds for quasilinear reaction-diffusion systems, Differ¬ 
ential Integral Equations, 8 (1995), 753-796. 

[53] Y. Tao, Global existence of classical solutions to a predator-prey model with non¬ 
linear prey-taxis, Nonlinear Anal. Real World Appl., 11 (2010), 2056-2064. 

[54] Y. Tao and Z.-A. Wang Competing effects of attraction vs. repulsion in chemotaxis. 
Math. Models Methods Appl. Sci., 23 (2013), 1-36. 

[55] T. Tona and N. Hieu, Dynamics of species in a model with two predators and one 
prey. Nonlinear Anal., 74 (2011), 4868-4881. 

[56] P. Turchin, Quantitative Analysis of Movement, Sinauer, Sunderland, Mass., 1998. 

[57] Q. Wang, C. Gai and J. Yan, Qualitative analysis of a Lotka-Volterra competition 
system with advection. Discrete Contin. Dyn. Syst., 35 (2015), 1239-1284. 

[58] [10.1007/s00332-016-9326-5] Q. Wang, Y. Song and L. Shao, Nonconstant positive 
steady states and pattern formation of ID prey-taxis systems, J. Nonlinear Sci., 
(2016), 1-27. 

[59] Q. Wang, J. Yang and L. Zhang, Time periodic and stable patterns of a 
two-competing-species keller-segel chemotaxis model: effect of cellular growth, 
preprint, http://arxiv.org/abs/! 505.06463. 



Patterns of systems with prey-taxis 


47 


[60] Q. Wang, L. Zhang, J. Yang and J. Hu, Global existence and steady states of a two 
competing species Keller-Segel chemotaxis model, Kinet. Relat. Models, 8 (2015), 
111-SOI. 

[61] X. Wang, W. Wang and G. Zhang, Global bifurcation of solutions for a predator- 
prey model with prey-taxis. Math. Methods Appl. Sci., 38 (2015), 431-443. 

[62] X. Wang and Q. Xu, Spiky and transition layer steady states of chemotaxis systems 
via globed bifurcation and Helly’s compactness theorem, J. Math. Biol., 66 (2013), 
1241-1266. 

[63] M. Winkler, Aggregation vs. global diffusive behavior in the higher-dimensioned 
Keller-Segel model, J. Differential Equations, 248 (2010), 2889-2905. 

[64] D. Xiao and S. Ruan, Codimension two bifurcations in a predator-prey system with 
group defense, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 11 (2001), 2123-2131. 

[65] J. Zhou, C.-G. Kim and J. Shi, Positive steady state solutions of a diffusive Leslie- 
Gower predator-prey model with Holling type II functioned response eind cross¬ 
diffusion, Discrete Contin. Dyn. Syst., 34 (2014), 3875-3899. 



